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Abstract 


This  research  postulates  and  demonstrates  a  modification  incorporating  rate- 
limited  sorption  effects  in  the  USGS  SUTRA  code  for  cleanup  of  a  hypothetical  sandy 
aquifer  by  pump-and-treat  remediation  methods.  Contaminant  transport  is  assumed  to  be 
affected  by  advection,  dispersion,  and  rate-limited  sorption/desorption.  Sorption  is 
assumed  to  be  either  equilibrium  or  rate-limited,  with  the  rate-limitation  described  by 
either  a  first-order  law,  or  by  Fickian  diffusion  of  contaminant  through  a  spherical 
immobile  pore  region.  Solutions  are  arrived  at  by  split  operator  methods  for  the  transport 
and  one-dimensional  Galerkin  solutions  for  the  solute  concentration  equations.  The 
resulting  model  is  tested  against  an  analytical  Laplace  transform  model  for  both  first-  and 
Fickian  diffusion  methods  in  a  radial  pumping  simulation.  Model  simulations  are  used  to 
evaluate  equilibrium,  first-order  and  Fickian  diffusion  effects  for  pulsed  and  continuous 
pumping  solutions  within  a  hypothetical  aquifer.  These  show  that  equilibrium  methods 
under-predicted  rebound,  while  first-order  methods  may  under-  and  over-predict  rebound 
within  the  matrix  for  certain  regions  and  may  be  equivalent  to  Fickian  diffusion  in 
equilibrium  regimes  for  cleanup  time  prediction.  Model  simulations  are  then  used  to 
show  the  efficiency  of  pulsed  pumping  methods  in  cleanup  mass  extraction  per  pumped 
volume  for  a  contaminated  aquifer  pump-and-treat  remediation  activity  versus  more 
conventional,  continuous  pumping  methods. 
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MODELING  OF  GROUND  WATER  AQUIFER  REMEDIATION 
BY  PULSED  PUMPING  WHEN  CONTAMINANT  TRANSPORT 

IS  AFFECTED  BY 

PHYSICAL,  NON-EQUILIBRIUM  SORPTION  AND  DESORPTION 

L  Introduction 

Background 

Ground  water  is  the  source  of  drinking  water  for  approximately  4b,000 
communities  and  twelve  million  individuals  across  the  country.  Almost  all  rural 
households  (95%)  and  50%  of  the  general  population  depend  upon  ground  water  as  their 
primary  drinking  water  source  [Jensen,  1993:250].  Spills  or  leaks  of  organic  chemicals  to 
the  environment  frequently  result  in  the  contamination  of  subsurface  soils  and  ground 
water  [Powers  and  others,  1991:463],  The  1989  Toxic  Release  Inventory  by  the 
Environmental  Protection  Agency  (EPA)  concluded  that  1 1%  of  toxic  chemical  releases 
occurred  on  land.  More  significantly,  14%  of  the  releases  were  underground  [U.S.  EPA, 
1989b:  173].  The  widespread  contamination  of  ground  water  due  to  years  of  accidental  or 
deliberate  disposal  of  various  chemicals,  particularly  where  contaminant  plumes  are 
encroaching  on  production  or  sole-source  aquifers  for  ground  water,  is  a  growing  issue 
within  the  United  States. 

Ground  water  contamination  exists  at  more  than  85%  of  the  1208  sites  included  on 
the  National  Priority  List  (NPL).  Further,  over  33,000  other  sites  have  been  identified  and 
included  in  the  Comprehensive  Environmental  Response,  Compensation,  and  Liability 
Information  System  for  ranking  and  potential  inclusion  on  the  NPL.  Contaminant  release 
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to  ground  water  has  been  identified  or  is  suspected  at  more  than  1700  Resource 
Conservation  and  Recovery  Act  facilities.  These  sites,  in  addition  to  the  large  number  of 
leaking  underground  storage  tank  sites  and  hazardous  waste  landfills  with  unlined  or 
ruptured  liners,  has  resulted  in  an  enormous  (ground  water)  contamination  problem  [Olsen 
and  Kavanaugh,  1993:42].  In  response,  the  Air  Force  is  engaging  in  a  program  to 
identify,  assess,  and  remediate  hazardous  waste  sites  at  military  installations  throughout 
the  United  States  [Goltz,  1991b:24]. 

Laboratory  researchers,  hydro-geologists  and  soil  scientists  have  actively  pursued 
modeling  methods  and  technologies  for  remediation  of  ground  water  for  years.  The  ability 
to  predict  rate  of  solute  migration  in  ground  water  systems  is  critical  to  the  evaluation  of 
potential  effects  on  water  quality  resulting  from  subsurface  disposal  of  hazardous  wastes, 
as  well  as  in  the  design  of  measures  to  restore  water  quality  of  contaminated  aquifers 
[Bahr  and  Rubin,  1987:438].  Due  to  the  rising  costs  in  performing  aquifer  remediation,  it 
is  important  that  decisions  made  regarding  ground  water  cleanup  be  based  upon  the  best 
available  information  [Adams  and  Virmontes,  1993:1-2].  One  major  source  of 
information  provided  to  the  Installation  Restoration  Program  (IRP)  decision  maker  comes 
from  contaminant  transport  models  [Goltz,  1991b:24]. 

Almost  all  remediation  of  ground  water  at  contaminated  sites  is  based  on  ground 
water  extraction  by  wells  or  drains  [Mackay  and  Cherry,  1989:630;  Goltz,  1 99  lb:25; 

Haley  and  others,  199 1 ;  1 19;  Hall  and  Johnson,  1992:216].  This  is  usually  accompanied  by 
treatment  of  the  extracted  water  prior  to  discharge.  Several  trends  in  the  effectiveness  of 
this  process  have  been  observed  in  the  field  over  time  [Mackay  and  Cherry,  1989:630; 
Haley  and  others,  1991:1 19;  Olsen  and  Kavanaugh,  1993:44]: 
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(1)  containment  of  ground  water  plumes  was  usually  achieved; 

(2)  contaminant  concentrations  dropped  significantly,  initially,  followed  by 
a  leveling  out  and; 

(3)  after  the  period  of  initial  rapid  decline,  the  continued  decreases  in 
concentration  were  usually  slower  than  anticipated  (sometimes  for  decades). 

Contaminant  transport  is  usually  modeled  by  the  advection/dispersion  equation 
[Goltz,  1986a:  1]  in  the  mobile  fluid  region,  with  a  diffusion  model  often  representing 
contaminant  transport  within  the  immobile  flow/soil  matrix  region.  The  transfer  of 
contaminant  between  the  flowing  ground  water  and  the  soil  matrix  is  modeled  by  the 
addition  of  a  sink  term  in  the  advection/dispersion  equation.  This  process  will  be  referred 
to  as  "sorption"  [Miller  and  Weber,  1988:466]. 

The  model  is  significantly  simplified  if  one  assumes  equilibrium  exists  between  the 
mobile  and  immobile  regions  within  the  soil  matrix.  Most  theoretical  and  mathematical 
models  make  an  assumption  that  local  equilibrium  (LEA)  is  occurring  within  the  macro  or 
micro-pores  of  the  soil  [Brusseau,  1992b:353].  In  this  case  the  solute  is  assumed  to  be 
sorbed  in  an  instantaneous,  linear  and  reversible  manner  [Goltz,  1986a:2].  If  the 
microscopic  processes  are  "fast  enough"  with  respect  to  the  bulk  flow  rate,  then  reversible 
sorption  reactions  are  assumed  to  be  in  a  state  of  local  equilibrium  [Valocchi,  1985:808], 

Models  used  to  predict  transport  of  solutes  that  assume  equilibrium  sorption  can 
predict  movement,  in  most  cases,  in  large  scale  aquifers  where  bulk  flow  rates  are  slow 
[Szcesody,  1988: 1].  However,  other  authors  have  demonstrated  that  a  LEA  model  may 
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significantly  underestimate  the  size  of  a  plume  of  reacting  solute  for  which  the  rate  of 
sorption  is  limited  by  non-equilibrium  mass  transfer  [Bahr  and  Rubin,  1987:438]. 

Several  laboratory  and  field  studies  have  indicated  that  LEA  may  not  be  adequate  for 
simulating  contaminant  transport  under  some  flow  conditions  [Valocchi,  1985:808;  Yu, 
1985:355;  Valocchi,  1986;1696,  Goltz  and  Roberts,  1986c:  1 139;  Mackay  and  others, 
1986:2047;  Brusseau  and  Rao,  1989d:41;  Brogan,  1991:165;  Goltz,  1991b:24;  Weber  and 
others,  1991:505],  While  local  equilibrium  may  exist  between  the  mobile  and  immobile 
zone  concentrations  in  soil  pores  under  natural  flow  fields,  the  increase  in  pore  water 
velocities  associated  with  pumping  of  the  aquifer  can  often  introduce  or  increase 
significantly,  non-equilibrium  conditions  [Brogan,  1991:2], 

Non-equilibrium  sorption  arises  when  physical  and  chemical  processes  at  the  single 
pore  level  are  slow  relative  to  advection  in  the  bulk  media  [Szecsody,  1988:1],  The 
chemical  binding  and  release  of  solute  onto  porous  and  solid  particles  within  the  soil 
matrix  cause  this  to  occur.  The  rate  of  sorption  is  controlled  by  one  of  three  possible  rate- 
limiting  steps  [Szecsody,  1988:9]: 

(1)  diffusion  through  immobile  fluid, 

(2)  diffusion  through  organic  matter,  and, 

(3)  chemical  binding  and  release. 

Non-equilibrium  models  have  been  introduced  to  better  characterize  these 
conditions.  These  models  are  based  on  the  assumption  that,  either  for  physical  or 
chemical  reasons,  adsorption  does  not  proceed  at  an  equal  rate  in  all  parts  of  the  soil 
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medium  [Van  Genuchten,  1981  :ii].  Conceptual  models  of  non-equilibrium  sorption  in 
soils  and  aquifers  are  grouped  into  two  classes:  (1)  physical  non-equilibrium  models, 
where  the  overall  sorption  rate  is  controlled  by  the  rate  at  which  the  solute  is  transported 
to  and  from  the  reacting  soil  surfaces  and  (2)  chemical  non-equilibrium  models  where  the 
overall  sorption  rate  is  equal  to  the  rate  of  reaction  at  the  soil-solution  interfaces 
[Valocchi,  1986:1694;  Grove  and  Stollenwerk,  1987:602-604]. 

The  effects  of  non-equilibrium  conditions,  'tailing'  and  'rebound',  have  been 
observed  in  several  laboratories  and  field  studies.  Tailing  describes  the  asymptotic 
decrease  in  the  rate  of  reduction  of  contaminant  concentration  in  extracted  water,  after  a 
relatively  rapid  initial  decrease.  Tailing  is  attributed  to  the  slow  diffusion  of  solute  into 
zones  of  immobile  water  [Van  Genuchten,  1981:3;  Yu,  1985:355;  Goltz,  1986a:5; 
Valocchi,  1986:1696;  Goltz  and  Roberts,  1988:37].  Rebound  involves  the  increase  in 
contaminant  concentration  observed  after  cessation  of  pumping  [Adams  and  Viramontes, 
1993:1-4].  Rebound  has  been  observed  several  years  after  the  pump-and-treat  system 
have  been  stopped  and  the  hazardous  site  closed  [Valocchi,  1986: 1694;  Mackay  and 
Cherry,  1989:633;  Goltz,  1991b:25].  Under  these  conditions,  it  appeared  that 
contaminant  in  the  sorbed  and  aqueous  phases  did  not  instantaneously  equilibrate,  but 
slowly  reached  equilibrium  [Goltz  and  Oxley,  1991a:547]. 

Non-equilibrium  solute  transport  is  recognized  as  a  key  limitation  to  contaminant 
removal  from  ground  water  by  "pump  and  treat"  methods  in  recent  literature  [Valocchi. 
1986:1699],  Slow  sorption  (with  attendant  tailing  and  rebound  effects)  of  organic 
substances  within  the  soil  has  resulted  in  long  periods  of  treatment,  sometimes  in  excess  of 
that  predicted  by  conventional  modeling,  with  attendant  costs.  Additionally,  such 
scenarios  lead  to  potentially  difficult  pumping  strategy  determinations.  In  many  situations, 
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continuous  pumping  did  not  prove  to  be  practical  [Mackay  and  Cherry,  1989:630;  Haley 
and  others,  1991: 1 19;  Olsen  and  Kavanaugh,  1993:44].  Processes  observed  to  improve 
the  plume  capture  and  reduce  quantities  of  water  extracted,  to  reach  the  goal  of  "clean 
water,"  more  efficiently,  are  referred  to  as  "pulse  pumping."  Field  experiments  and 
modeling  studies  have  demonstrated  that  this  optimization  of  the  process  may  achieve 
favorable  results  [Brogan,  1991:165;  Haley  and  others,  1991:122;  Goltz,  1991b:25; 
Haggerty,  1992:37;  Adams  and  Virmontes,  1993:5-3]. 

An  extensive  survey  [Van  der  Heijde  and  Elnawawy,  1993]  of  142-mathematical 
models  for  ground  water  transport  showed  that  51  models  use  some  form  of  an 
equilibrium  or  linear  sorption  routine  to  compute  solute  movement,  and  only  three 
evaluate  non-equilibrium  sorption.  Huso  (1989),  using  finite  elements,  and  Adams  and 
Virmontes  (1993)  using  an  analytical  model,  modeled  physical,  non-equilibrium  in  radial, 
pulsed  pumping  tests.  Brogan  [1991:27]  extensively  modeled  non-equilibrium  solute 
transport  and  pulsed  pumping  for  chemical  (kinetic)  non-equilibrium  using  a  two- 
dimensional,  reversible,  first-order  kinetic  solute  transport  module  addition  to  the  public 
domain  USGS  code,  SUTRA.  Brogan's  work  was  extended  for  multi-contaminant 
transport  using  modified,  first-order  kinetic  computation  equations,  for  faster  computation 
[Haggerty,  1992:42-44], 

SUTRA  is  a  two-dimensional  public  domain  code,  which  calculates  subsurface 
ground  water  movement  and  contaminant  transport  in  both  the  saturated  and  unsaturated 
zones.  It  models  single  species  transport  through  a  quadrilateral  hybrid  finite  element 
transport  model.  Spatial  derivatives  are  approximated  with  finite  elements  and  temporal 
derivatives  with  finite  differences  [Voss,  1984:5-8].  A  modification  to  SUTRA,  for 
physical  non-equilibrium,  utilizing  the  two-region  model  equations  developed  by  Van 
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Genuchten  [1985:513]  for  both  first-  and  second-order,  non-equilibrium  sorption  will  be 
investigated  as  a  more  rigorous  approach  towards  modeling  the  physical  processes  within 
the  soil  than  the  chemical  approach  presented  by  Brogan. 

Specific  Problem 

The  purpose  of  this  research  is  to  mathematically  model  aquifer  remediation  by 
pulsed  pumping  when  contaminant  transport  is  affected  by  physical,  diffusion,  rate-limited 
sorption.  This  research  will  extend  the  work  of  Goltz  and  Oxley  (1991a),  Brogan  (1991) 
and  Adams  and  Virmontes  (1993). 

Research  Objectives 

The  specific  objectives  of  this  research  are  to: 

1.  Develop  and  validate  a  modification  to  the  SUTRA  transport  model  that  incorporates 
the  mechr  ..isms  of  advection,  dispersion,  and  physical  non-equilibrium  sorption  using  a 
two-region  model  to  emulate  immobile  regions  within  the  soil  matrix  by  multiple 
simulations.  The  study  is  limited  to  the  special  case  where  advection  is  due  to  steady 
ground  water  velocity. 

2.  Perform  a  comparison  test  with  an  existing  analytical  code  for  both  first-  and  second- 
order,  rate-limited  sorption  for  a  radial  flow  zone,  simulated  by  a  pumping  well. 

3.  Compare  and  contrast  solute  extraction  efficiency  of  the  second-order,  physical  non¬ 
equilibrium  sorption  algorithm  with  first-order  and  an  LEA  model,  for  a  pulsed-pumping 
remediation  simulation. 
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4.  Compare  efficiency  of  contaminant  extraction  between  continuous  pumping  strategies 
and  pulsed  pumping  strategies  through  integration  of  breakthrough  curves  during 
extraction  and  comparison  of  contaminant  mass  extracted  with  the  pumped  water. 

Scope  and  Limitations 

Simplifying  assumptions  used  in  this  research  are  listed  below: 

1.  Single  species,  solute  transport  will  be  simulated. 

2.  Boundary  conditions  are  constant  throughout  pumping  and  rest  cycles. 

3.  Transport  of  the  contaminant  due  to  the  natural  ground  water  gradient  is  negligible 
during  the  period  of  extraction.  Containment  of  the  plume  is  not  necessary  to  prevent  the 
spread  of  ground  water  contamination  to  other  areas  (outside  the  representative  volume  to 
be  studied).  The  extraction  well  must  create  a  flow  field  that  encompasses  the  area  of  the 
contaminant  plume. 

4.  In  the  governing  contaminant  transport  equation,  molecular  diffusion  is  negligible  with 
respect  to  mechanical  dispersion,  when  the  pump  is  on. 

5.  The  formation  is  homogenous.  The  vertical  dimension  is  a  fixed  parameter.  This 
facilitates  SUTRA's  three-dimensional  characterization  of  the  aquifer  and  plume. 

6.  The  maximum  concentration  is  located  at  the  center  of  the  plume,  with  the  radius  of 
the  plume  defined  as  the  distance  at  which  the  concentration  is  1%  of  the  minimum  action 
level. 
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7.  Hydraulic  head  and  contaminant  concentration  are  both  constant  with  respect  to  depth. 


8.  All  water  entering  the  aquifer  from  the  boundaries  is  contaminant  free. 

9.  It  is  not  desirable  to  de-water  the  aquifer,  nor  to  pump  the  well  dry. 

Definitions 

Key  terms  associated  with  contaminant  transport  and  aquifer  remediation,  as 
defined  by  the  EPA  unless  otherwise  specified,  are  listed  below  [Mercer  and  others,  1990; 
U.S.  EPA,  1993]. 

1.  Absorption:  adherence  of  ions  or  molecules  in  solution  to  the  surface  of  solids. 

2.  Advection:  the  process  whereby  solute  is  transported  by  the  bulk  mass  of  flowing 
fluid. 

3.  Aquifer:  a  geologic  unit  that  contains  sufficient  saturated  permeable  material  to 
transmit  significant  quantities  of  water. 

4.  Aquitard:  a  relatively  impermeable  layer  that  greatly  restricts  the  movement  of  ground 
water  [Freeze  and  Cherry,  1979:47]. 

5.  Breakthrough  Curve  (BTCl:  plot  of  contaminant  concentration  versus  time  relation  at 
a  particular  point  (spatially)[Freeze  and  Cherry,  1979:391], 
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6.  Cleanup:  the  attainment  of  a  specified  contaminant  concentration  level  [Goltz  and 
Oxley,  1991a:547];  water  quality  or  Safe  Drinking  Water  Act  Standard  is  sometimes  used 
as  the  regulatory  criteria  to  define  cleanup  of  an  aquifer.  This  concentration  level  is  used 
to  define  the  solute  concentration  criteria  to  define  cleanup  for  this  thesis.  It  is  also 
necessary  for  solute  concentrations  to  remain  below  this  cleanup  criterion  after  pumping 
has  stopped. 

7.  Cleanup-mass:  total  mass  of  contaminant  extracted  by  a  pump-and-treat  system  as 
dissolved  solute  within  the  pumped  ground  water. 

8.  Cleanup-time:  total  time  it  takes  to  achieve  cleanup  [Brogan,  1991:82]. 

9.  Cleanup- volume:  total  volume  of  ground  water  that  is  extracted  from  the  aquifer 
[Brogan,  1991:82]. 

10.  Concentration  Gradient:  movement  of  a  contaminant  from  a  region  of  higher 
concentration  to  one  of  lower  concentration  [Freeze  and  Cherry,  1979:25]. 

11.  Contaminant  Removal  Efficiency:  the  amount  of  contaminant  removed  per  volume  of 
water  pumped.  This  is  the  slope  of  the  mass  arrival  curve  [Brogan,  1991:40]. 

12.  Desorption:  the  reverse  of  sorption. 

13.  Diffusion:  mass  transfer  as  a  result  of  random  motion  of  molecules.  It  is  described  by 
either  a  first-order  equation  of  Ficks  law  of  diffusion  (second-order  equation). 
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14.  Dispersion:  the  spreading  and  mixing  of  the  contaminant  in  ground  water  caused  by 
diffusion  and  mixing  due  to  microscopic  variations  in  velocities  within  and  between  pores. 

15.  Effective  Porosity:  the  ratio,  usually  expressed  as  a  percentage,  of  the  total  volume  of 
voids  available  for  fluid  transmission  to  the  total  volume  of  the  porous  medium. 

16.  Extraction  Well:  a  pumped  well  used  to  remove  contaminated  ground  water; 
operated  solely  during  the  pumping  cycle  and  has  the  primary  objective  of  removing 
contamination  from  the  aquifer.  The  extraction  well  provides  plume  containment  during 
the  pumping  cycle  [Brogan,  1991: 144]. 

17.  Gradient  Control  Well:  a  pumped  well  used  to  remove  contaminated  ground  water 
during  the  rest  cycle  with  the  objective  of  preventing  off-site  migration  of  the  contaminant 
plume  while  desorption  takes  place  [Brogan,  1991:144]. 

18.  Heterogeneous:  A  geologic  unit  in  which  the  hydrologic  properties  vary  from  point 
to  point. 

19.  Homogeneous:  a  geologic  unit  in  which  the  hydrologic  properties  are  identical  from 
point  to  point. 

20.  Mass  Arrival  Curve:  a  plot  of  the  mass  fraction  removed  versus  the  number  of  pore 
volumes  pumped  [Brogan,  1991:40]. 

2 1 .  Pulsed  Pumping:  a  pump-and-treat  enhancement  where  extraction  wells  are 
periodically  not  pumped  to  allow  concentrations  in  the  extracted  water  to  increase. 
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22.  Pumping  Cycle:  time  period  during  which  the  extraction  well  is  operating  to  remove 
contaminant  mass  from  the  ground  water  [Brogan,  1991 : 108]. 

23.  Rebound:  increase  in  contaminant  concentration  in  the  ground  water  that  is  observed 
due  to  desorption  [Adams  and  Virmontes,  1993:1-4]. 

24.  Rest  Cycle:  non-pumping  period  for  the  extraction  well.  It  may  include  a  pumping 
period  for  a  gradient  control  well,  where  the  objective  is  to  allow  solute  concentrations  to 
increase  through  desorption  [Brogan,  1991:  108]. 

25.  Retardation:  the  movement  of  a  solute  through  a  geologic  medium  at  a  velocity  less 
than  that  of  the  flowing  ground  waters  due  to  sorption  or  other  removal  of  the  solute. 

26.  Sorption:  processes  that  remove  solutes  from  the  fluid  phase  and  concentrate  them 
on  the  solid  phase  of  a  medium;  the  generic  term  used  to  encompass  the  phenomena  of 
adsorption  and  adsorption. 

27.  Tailing:  the  slow,  nearly  asymptotic  decrease  in  contaminant  concentration  in  water 
flushed  through  contaminated  material. 

28.  Unconfined  Aquifer:  an  aquifer  in  which  the  water  table  forms  the  upper  boundary 
[Freeze  and  Cherry,  1979:48]. 
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Overview 

IRP  planners  use  contaminant  transport  models  to  assess  risk;  to  design  remedies 
and  to  estimate  remediation  cost  and  cleanup  time  duration  at  IRP  hazardous  sites. 

Chapter  I  examined  one  modeling  assumption  often  employed  at  these  sites.  It  discussed 
how  this  assumption  does  not  account  for  slow  or  non-equilibrium  sorption/desorption 
that  has  been  observed  in  several  field  and  laboratory  studies.  A  modification  to  the 
SUTRA  program  was  proposed  to  evaluate  a  pulsed-pumping  scheme  as  a  technique  to 
enhance  the  effectiveness  of  ground  water  pump-and-treat  systems  by  accounting  for 
physical  non-equilibrium  sorption  in  the  soil  matrix.  This  chapter  concludes  with  a 
research  proposal  to  validate  a  mathematical  model  of  an  aquifer  decontamination  by  using 
pulsed-pumping  when  the  contaminant  transport  is  affected  by  physical  non-equilibrium 
sorption  and  desorption. 

Chapter  II  discusses  the  literature  associated  with  sorbing  solute  transport 
modeling.  An  introduction  to  the  basic  equation  used  to  represent  solute  transport  in 
ground  water  is  presented.  Then,  the  chapter  reviews  the  efforts  of  researchers  to  develop 
mathematical  models  to  account  for  equilibrium  and  non-equilibrium  sorption.  A 
discussion  follows  on  the  validity  of  the  respective  regimes  (equilibrium,  versus  non¬ 
equilibrium).  Physical  versus  chemical  (kinetic)  non-equilibrium  modeling  of  processes  is 
then  discussed.  Finally,  a  discussion  concerning  the  application  of  various  pump-and-treat 
technologies  (continuous  pumping  versus  pulsed-pumping)  is  presented. 

Chapter  III  presents  a  derivation  of  the  Van  Genuchten  [1985]  two-region  model, 
as  modified  by  Goltz  and  Oxley  [Goltz  and  Oxley,  1991a;  Goltz  and  Roberts,  1986b,  c] 
for  physical  non-equilibrium  and  its  application  in  modification  of  the  SUTRA  code  to 
accommodate  first-  and  second-order,  rate-limited  sorption  within  spherical  aggregates. 
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These  describe  the  contaminant  transport  by  means  of  radial  flow  with  sorption  in  an 
aquifer  undergoing  pulsed-pumping  operations.  Sorption/desorption  is  described 
assuming  Fickian  diffusion  of  the  contaminant  through  immobile  water  in  zones  of 
spherical  geometry. 

In  Chapter  IV,  we  validate  the  modifications  made  to  the  SUTRA  code  for  finite 
elements  and  differences  with  a  previously  validated  analytical  model  [Adam's  and 
Viramontes,  1993]  which  incorporates  rate  limited  sorption  and  desorption.  Breakthrough 
times,  approach  to  equilibrium,  rebound  and  cleanup  time  are  compared  for  first-  and 
second-order  simulations  with  an  LEA  model  for  a  range  of  parameters.  Several 
breakthrough  curves  are  generated  and  used  to  illustrate  the  simulations  and  model 
comparisons.  We  then  analyze  simulations  of  pulsed  and  continuous  pumping  schemes  for 
optimal  extraction  efficiency.  Cleanup  time,  mass  and  pore  volumes  extracted  are 
computed  and  compared  for  pulsed,  pump-and-treat  remediation  activities  versus 
continuous  pumping  strategies  in  remediation  of  a  hypothetical  sandy  aquifer. 

Chapter  V  summarizes  the  research,  draws  conclusions  based  on  the  findings  and 
offers  recommendations  for  further  research. 
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II.  Literature  Review 


Overview 

Hazardous  waste  remediation  will  become  even  more  important  in  the  future  than 
at  the  present.  As  of  1970,  the  EPA  estimated  that  94%  of  the  rural  population  and  37% 
of  the  public  water  supplies  relied  on  ground  water  as  their  primary  source  [Canter  and 
others,  1987:7].  An  Office  of  Technology  Assessment  (OTA)  estimate  of  major  sources 
of  potential  ground  water  contamination  include:  hazardous  waste  sites  (7,000);  leaking 
underground  storage  tanks  (75,000-100,000);  industrial  and  municipal  landfills  (100,000) 
and  septic  tank  systems  (20  million)  [Canter  and  others,  1987:67-90]. 

The  United  States  is  currently  spending  approximately  $10  billion  annually  on 
environmental  cleanup  and  this  is  expected  to  rise  to  $1 1.2  billion  by  1999  for  the  Federal 
sector  alone  [Bredehoeft,  1994:97].  The  1992  Department  of  Defense  (DOD) 
Authorization  Act  provided  over  $1.4  billion  for  environmental  restoration  for  all  services. 
The  Air  Force  has  estimated  that  $1.45  billion  will  be  required  to  implement  existing 
environmental  interim  agreements  for  remediation  projects  currently  identified  [DOD 
DERP,  1992:27,  A-5], 

A  recent  study  of  the  Superfund  National  Priorities  List  (NPL)  found  that  46%  of 
the  NPL  sites  used  pump-and-treat  technology  for  collection  and  treatment  [Bredehoeft, 
1994:98]  and  13%  of  DOD  Installation  Restoration  Program  (IRP)  sites  involved  ground 
water  treatment  [DOD  DERP,  1992:8].  The  Air  Force  operates  331  installations  and 
owns  13  contractor-operated  plants,  within  which  there  are  3520  IRP  sites  currently  open 
for  remedial  investigation  or  action  [DOD  DERP,  1992:6].  As  outlined  above,  many  of 
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the  problems  endemic  to  the  civilian  sector  are  equally  prevalent  at  Air  Force  installations. 
The  four  major  source  types  of  pollution  of  ground  water  for  Air  Force  installations  are 
[Gorelick  and  others,  1993:32]: 

(1)  leaking  underground  storage  tanks; 

(2)  leaking  pipelines  (including  sewers); 

(3)  spills  or  other  major  releases  of  liquids  at  the  ground  surface; 

(4)  industrial  operations. 

The  most  common  contaminants  at  Air  Force  sites  which  are  being  cleaned  up 
under  the  IRP  comprise  three  major  contaminant  groups  [Gorelick  and  others,  1993:2-3]: 

(1)  petroleum  hydrocarbon  contaminants,  which  usually  include  benzene,  ethyl 
benzene,  toluene  and  xylene, 

(2)  trichloroethylene  (TCE)  and  related  solvents  from  degreasing,  paint  handling, 
metal  plating  and  other  similar  operations, 

(3)  toxic  metals  and  other  inorganic  contaminants  such  as  lead,  cadmium, 
chromium  and  cyanide. 

Given  the  ground  water  pollution  problems  in  the  civilian  sector  and  similar  ones 
being  experienced  at  DOD  installations  under  the  IRP,  it  is  important  that  modeling  efforts 
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continue  to  better  quantify  and  clarify  the  spatial  and  temporal  characteristics  of 
contaminant  plumes,  prior  to  initiation  of  expensive  remediation  projects.  Accurate 
modeling  of  sorption  processes  is  essential  in  assessing  the  effects  of  transport  of  sorbing 
organic  solutes  for  optimization  of  pump-and-treat,  treatment  technology.  According  to 
recent  estimates  for  the  Superfund  Program,  the  average  per-site  cost  of  characterization 
through  initiation  of  remedial  action  for  a  pump-and-treat  system  is  about  $14.1  million 
[Lee  and  Kitanidis,  1991:2203].  Given  the  expense  of  remediation,  risk  to  public  health, 
and  time  involved  for  clean-closing  of  a  pump-and-treat  site;  the  ability  to  model 
contaminant  transport,  both  spatially  and  temporally  may  give  critical  insight  towards 
design  strategies  which  optimize  the  process  of  remediation. 

Critical  to  this  is  the  development  of  a  mathematical  model  which  can  accurately 
simulate  solute  transport  through  the  subsurface.  Of  the  various  transport  processes, 
sorption  is  one  of  the  most  important;  it  is  often  the  paramount  process  controlling  the 
behavior  of  contaminants  in  the  sub-surface  under  many  conditions  [Piwoni  and  Keeley, 
1990: 1].  Sorption  results  in  solutes  traveling  slower  than  the  water  present  in  the  aquifer 
[Thorbjamarson  and  Mackay,  1994:401].  If  sorption  is  linear,  reversible  and 
instantaneous,  the  solute  transport  can  be  modeled  using  simplistic  methods,  such  as 
equilibrium  transport.  These  simplifying  assumptions  are  called  the  "Local  Equilibrium 
Assumption  (LEA)".  However,  these  assumptions  are  not  justified  under  all 
circumstances  and  non-equilibrium  transport  of  organic  solutes  during  pump-and-treat 
remediation  of  contaminated  aquifers  has  lead  to  tailing  of  concentrations  above  typical 
clean-up  goals  for  very  long  periods  of  time  [Mackay  and  Cherry,  1989:630], 

The  effect  of  contaminant  desorption  from  solids  is  to  extend  the  time  required  to 
remove  specified  levels  of  contaminant  for  a  given  pumping  rate  due  to  retardation  of  the 
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contaminant  plumes  migration.  Further,  this  retardation  may  require  extraction  of  greater 
quantities  of  water  than  estimated  under  the  LEA  case,  due  to  the  desorption  of 
contaminant  from  immobile  zones  (through  diffusion  processes)  to  mobile  zones  where 
"clean  water"  has  replaced  that  removed  during  pumping  [Mackay  and  Cherry,  1989:633; 
Hall  and  Johnson,  1992:218]. 

The  analytical  and  mathematical  modeling  of  the  processes  for  non-equilibrium 
transport  of  organic  solutes  has  received  considerable  attention  over  the  past  five-years. 
Modeling  the  processes  that  lead  to  non-equilibrium  transport  of  organic  solutes  under 
field  conditions,  particularly  forced-gradient  pumping  conditions  such  as  those  induced  in 
most  remediation  efforts  has  received  some  emphasis  of  late  from  some  researchers, 
notably  Huso  (1989),  Brogan  (1991),  Adam's  and  Viramontes  (1993)  for  single  species 
transport  and  Haggerty  (1992)  for  multi-species  transport. 

Contaminant  Transport  Processes 

Sorption  Isotherm  Linearity 

Mass  transport  is  controlled  by  a  variety  of  physical,  chemical  and  biological 
processes.  Quantitative  descriptions  of  the  advection  and  dispersion  processes  are  well 
understood.  Research  emphasis  has  focused  on  diffusion  and  sorption  processes, 
attempting  to  model  these  processes  at  the  micropore  level,  and  to  account  for  the  effects 
of  tailing  (slow  desorption)  of  solute  over  long  periods  of  time  for  certain  solute  and  soil 
matrixes.  Traditionally,  three  approaches  have  been  taken  in  sorption  modeling  [Parker 
and  Valocchi,  1986:399]: 

(1)  bicontinuum  approach  with  interchange  between  mobile  and  immobile  zones  of 
specified  geometry  governed  by  Ficks  Law  for  diffusion, 
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(2)  bicontinuum  approach  with  interchange  described  by  an  empirical  first-order 
kinetic  expression, 

(3)  equivalent  monocontinuum  approach  which  assumes  residence  times  are 
sufficiently  large  to  achieve  negligible  concentration  differences  between  mobile  and 
immobile  pore  regions  (LEA). 

Figure  2. 1  illustrates  the  processes  which  are  assumed  to  occur  within  the 
micropore  zone  for  contaminant  transport: 


■  Tortuosity 

Diffusion  through  Bound 
Organic  Layer 


Figure  2.1.  Bicontinuum  Model  for  Contaminant  Transport 


The  typical  approach  in  modeling  contaminant  transport  has  been  to  assume 
conditions  of  local  equilibrium  [Bear,  1987:171;  Brusseau  and  Rao,  1989d:34],  The 
conventional  equation  used  to  describe  the  advective-dispersive,  one-dimensional 
transport  of  a  sorbing  solute  is  given  by  [Lapidus  and  Amundson,  1952:985;  Van 
Genuchten,  1981:2]: 
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(2.1) 


ac  pas  ^a2c  ac 
at  0  at  ax2  ax 


where  C  is  the  solute  concentration  in  aqueous  phase,  0  is  the  volumetric  moisture 


content,  x  is  the  distance  in  flow  direction,  D  is  the  dispersion  coefficient,  p  is  the  bulk 


density  of  soil ,  S  is  the  mass  of  solute  adsorbed  per  unit  dry  mass  of  soil,  t  is  time  and 


v  is  the  average  fluid  velocity. 


Describing  the  dS/dt  term  is  often  the  critical  step  in  model  formulation.  The 
LEA  results  in  the  simplified  expression: 


as  ac 

Yt~Kd~dt 


(2.2) 


where  Kj  is  the  partition  coefficient  between  soil  and  water  [Szecsody,  1988:3]. 
Upon  substitution,  the  transport  equation  becomes: 


a2c 


R^-  =  D 
dt  dx 


ac 


2  v'a* 


(2.3) 


where 


/?  =  !  + 


0 


(2.4) 


R  is  the  retardation  factor.  Equation  (2.3)  will  hereafter  be  referred  to  as  the  advective- 
dispersive  (A-D)  transport  equation. 
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Advection  is  the  movement  of  dissolved  solute  with  flowing  ground  water.  It  is 
9C 

represented  by  the  v —  term  in  the  A-D  equation.  In  coarse-grained,  homogeneous 
dx 

aquifers,  advective  transport  dominates  the  transport  process.  In  the  design  of  capture 
and  containment  systems,  it  is  common  to  treat  advection  as  the  sole  mechanism  for 
contaminant  transport  [Gorelick  and  others,  1993:42]. 

The  two  processes  that  constitute  hydrodynamic  dispersion  are  molecular  diffusion 

and  mechanical  dispersion.  Molecular  diffusion  is  the  process  in  which  dissolved 

contaminants  move  from  areas  of  high  concentration  to  areas  of  low  concentration  in 

response  to  the  presence  of  a  concentration  gradient.  Mechanical  dispersion  is  the 

component  of  transport  that  results  from  the  convoluted  paths  that  water  and  contaminant 

particles  follow  while  flowing  through  porous  and  fractured  media  (reflected  by  a 

tortuosity  term) [Gorelick  and  others,  1993:42;  Freeze  and  Cherry,  1979:75).  Diffusion 

d^C 

and  dispersion  are  included  in  the  D- —  term  in  the  A-D  equation. 

dx2 

The  physical  significance  of  the  retardation  factor  (Equation  [2.4])  is  that  it 
measures  how  much  slower  the  solute  migrates  than  water  because  the  solute  spends  part 
of  its  time  sorbed  on  the  soil  matrix  and  is  immobile  [Charbeneau  and  others,  1992:43], 
The  higher  the  fraction  of  contaminant  sorbed,  the  more  retarded  is  its  transport.  If  the 
sorptive  reaction  behaves  linearly  and  is  at  equilibrium,  the  solute  will  move  at  an  average 
velocity  equal  to  the  average  linear  velocity  of  the  ground  water  flow  divided  by  the 
retardation  factor  [Gorelick  and  others,  1993:44]. 
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When  LEA  is  applied,  equilibrium  is  assumed  to  exist  at  each  point  in  space 
between  dissolved  solute  (associated  with  the  liquid  phase)  and  solute  (associated  through 
sorption),  with  the  solid  phases  of  the  medium  [Harmon  and  others,  1989:406],  and  is 
modeled  using  a  partition  coefficient  as  in  Equation  (2.2).  As  mentioned  previously,  this 
assumption  greatly  simplifies  the  transport  model. 

Linear  isotherms  are  commonly  used  to  model  sorption  in  solute  transport  at  low 
concentrations  by  assuming  that  a  linear  relationship  exists  between  dissolved  and  sorbed 
solute  [Harmon  and  others,  1989:407].  However,  in  many  cases,  linear  models  do  not 
adequately  describe  the  solute  relationship.  As  a  result,  several  non-linear  relationships 
have  been  applied  under  the  LEA.  The  commonly  used  isotherms  are  summarized  in 
Table  2.1  below: 


TABLE  2.1. 

ADSORPTION  ISOTHERMS 


Adsorption  Isotherm 


Equation 


Linear 

Freundlich 

Langmuir 

Langmuir  Two-Surface 
Competitive  Langmuir 


S  =  KrfC 
S  =  KCn 

S  =  kbC/(  1  +  kC) 

S  =  (kibiC/(l  +k1C))  +  (k2b2C/(l 
+k2C)) 

(C]/C2)/S  =  (b2/k|bj)  +  (C|/b2C2) 

fEngesgaard  and  Christensen.  1988:1921 


where  S  is  the  sorbed  species;  K<j  is  the  distribution  coefficient;  Kj  is  the  bond 


.strength;  Cj  is  the  dissolved  species  concentration;  kj  and  k2  are  the  forward  and 
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backward  reaction  rate  constants  [Carnahan  and  Remer,  1984:241]  and  bj  is  the 
maximum  capacity  of  sorbed  solute. 

For  many  nonionic  compounds  the  Freundlich  adsorption  isotherm  is  typically 
used  [Jury  and  Ghodrati,  1989:275;  Pignatello,  1989:46].  Non-linearity  is  important  for 
ionic  or  strongly  polar  compounds  and  is  commonly  modeled  using  the  Langmuir  isotherm 
[Brusseau  and  Rao,  1989d:38].  Linear  sorption  isotherms  are  frequently  used  where  the 
distribution  coefficient  for  the  chemical  species  has  been  found  to  be  a  function  of  the 
hydrophobic  character  of  the  organic  compound  and  the  amount  of  organic  matter 
present,  using  the  equation  [Charbeneau  and  others,  1992:42;  Szecsody,  1989:3;  Johnson 
and  others,  1989:43]: 


K,  =  (2.5) 

where  KqC  is  the  organic  carbon  partition  coefficient  and  ^oc  is  the  fraction  of  organic 
carbon  within  the  soil  matrix. 

The  tailing  effect  on  breakthrough  curves  results  in  a  "spreading  front"  due  to  the 
slower,  low  concentration  components  of  the  plume  lagging  behind  the  faster  moving, 
high  concentration  components  of  the  contaminant  plume.  The  tailing  effect  has  been 
further  explained  by  other  authors  [Hall  and  Johnson,  1992:218]  as  representing 
heterogeneities  in  permeability  within  the  soil  matrix,  where  pumping  initially  extracts 
contaminants  from  zones  where  advection  dominates  (larger  open  pores  where  fluid 
moves  more  readily)  and  where  solute  is  retained  in  the  smaller  closed  pores.  Over  time 
the  contaminant  in  the  closed  pores  is  postulated  to  migrate  to  equilibrium  with  the 
concentration  gradient  in  the  open  pores  by  a  process  of  diffusion.  This  recharge  of 
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contaminant  within  the  open  pores  then  causes  the  tailing  effect  explained  mathematically 
above.  The  tailing  effect  serves  to  prolong  the  expected  time  of  extraction  of  contaminant 
and  increases  the  volume  of  water  necessary  for  extraction  and  treatment  [Mercer  and 
others,  1990:8]. 


LEA  Validity 

As  previously  mentioned,  LEA  is  a  major  component  of  the  conventional  solute- 
transport  equation.  If  non-equilibrium  conditions  are  present,  LEA-based  models  will 
predict  a  breakthrough  response  that  occurs  too  late  and  exhibits  too  little  dispersion 
[Valocchi,  1985:808],  For  the  LEA  to  be  valid,  the  rate  of  the  sorption  process  must  be 
fast  relative  to  the  other  processes  affecting  solute  concentration  (e.g.,  advection, 
hydrodynamic  dispersion)  so  that  equilibrium  may  be  established  between  the  sorbing 
contaminant  and  the  pore  fluid  [Bahr  and  Rubin,  1987:438;  Engesgaard  and  Christensen, 
1988: 188;  Brusseau  and  Rao,  1989d:41].  Additionally,  other  authors  have  described 
when  the  LEA  applies  as: 

(1)  when  coefficient  of  hydrodynamic  dispersion  (D)  was  of  same  order  of 
magnitude  as  the  effective  molecular  diffusion  coefficient  (De)  [Valocchi,  1985:808;  Bahr 
and  Rubin,  1987:438]. 

(2)  when  the  limiting  (smallest)  Damkohler  I  (Da  I)  number  (an  empirical 
relationship  between  the  approach  to  local  equilibrium  and  values  of  dimensionless 
parameters  representing  ratios  of  an  average  transit  time  to  a  reaction  time)  exceeds  100 
[Bahr  and  Rubin,  1987:439;  Brusseau  and  others,  1991b:  1 140-1 141]. 
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(3)  when  the  chemical  reaction  of  the  solute  is  fast  with  respect  to  the  bulk  flow  of 
fluid  [Valocchi,  1985:808;  Grove  and  Stollenwerk,  1987:603;  Bahr  and  Rubin,  1987:438; 
Engesgaard  and  Christensen,  1988:188]. 

Sorption  Non-Equilibrium 

Sorption  data  from  batch  experiments  have  been  found  to  exhibit  a  two-stage 
approach  to  equilibrium:  a  short  initial  phase  of  fast  uptake,  followed  by  an  extended 
period  of  much  slower  uptake  [Wu  and  Gschwend,  1986:724].  This  was  also  observed 
during  the  Borden  field  test  where  long  term  increases  in  retardation  were  attributed  to 
non-linear  equilibrium  behavior  or  non-equilibrium  sorption-desorption  [Roberts  and 
others,  1986:2056;  Bahr  and  Rubin,  1987:438]. 

Data  from  many  laboratory  experiments  exhibit  asymmetrical  BTCs,  with  earlier 
breakthrough,  increased  time  to  reach  complete  breakthrough  (i.e.,  breakthrough-front 
tailing)  and  increased  time  for  complete  desorption  (i.e.,  elution-front  tailing).  A  rate- 
limited  non-equilibrium  process  has  been  postulated  as  being  responsible  for  the  early 
breakthrough  and  tailing  that  produces  the  asymmetry  exhibited  by  BTCs  influenced  by 
non-equilibrium  [Brusseau  and  Rao,  1989d:43J.  The  occurrence  of  sorption  non¬ 
equilibrium  increases  solute  dispersion.  Hydrodynamic  dispersion  will  tend  to  dominate 
under  conditions  of  low  pore-water  velocities  and  in  fine-grained,  poorly  aggregated 
media,  whereas  non-equilibrium  induced  dispersion  tend  to  predominate  for  larger 
velocities  and  for  large-scale  aggregation  [Brusseau  and  Rao,  1989d:43;  Rao  and  others. 
1980:1145]. 
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Non-Equilibrium  Sorption  Models 

Data  exhibiting  behavior  that  deviates  from  that  predicted  by  the  simple  LEA 
model  has  been  reported,  thus  challenging  the  validity  of  this  approach.  As  previously 
stated,  there  are  cases  where  equilibrium  is  not  reached  rapidly  enough  with  respect  to 
advective  transport  to  allow  use  of  A-D  equations  with  equilibrium  sorption.  As  outlined 
in  Chapter  I,  non-equilibrium  sorption  and  desorption  in  sub-surface  media  arises  when 
physical  and  chemical  processes  at  the  single  pore  level  are  slow  relative  to  advection  in 
the  bulk  media  [Szecsody,  1988: 1].  Non-equilibrium  sorption  has  been  attributed  to 
several  factors  [Brusseau  and  Rao,  1989d:33]: 

(1)  kinetic  sorption; 

(2)  diffusive  mass  transfer  resistance; 

(3)  sorption-desorption  non-singularity. 

Physical  Non-Equilibrium  Models 

Physical  non-equilibrium  is  defined  as  the  rate-limited  sorption/desorption  resulting 
from  the  existence  of  mass-transfer  resistance  [Brusseau  and  Rao,  1989d:62].  Some 
authors  have  postulated  physical  non-equilibrium  models  which  hypothesize  the  existence 
of  zones  of  immobile  water  which  act  as  diffusion  sources  and  sinks  for  solutes  [Van 
Genuchten,  1985:513;  Goltz  and  Roberts,  1986b:77;  Szecsody,  1988:8;  Brusseau  and 
Rao,  1989d:45].  Further,  these  models  can  be  described  as  either  first-order  rate  models 
[Van  Genuchten  and  Wierenga,  1976:474;  Goltz  and  Roberts,  1986b:81]  or  second-order 
diffusion  models  [Van  Genuchten  and  Wierenga,  1976:474;  Rao  and  others,  1980: 1 139; 
Van  Genuchten  and  others,  1984:335;  Brusseau  and  Rao,  1989d:47].  For  hydrophobic. 
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organic  compounds,  physical  non-equilibrium  and  intra-organic  matter  diffusion  have  been 
identified  as  probable  predominant  factors  causing  non-ideality  within  the  organic  matter 
matrix  [Brusseau,  1989a:vi]. 

Most  physical  non-equilibrium  models  incorporate  a  two-zone  description  of  a 
porous  medium,  although  some  authors  have  postulated  more  regions.  This  divides  the 
medium  into  an  immobile  zone  and  a  mobile  zone  into  which  the  solute  diffuses. 

Transport  in  the  mobile  zone  is  governed  by  the  advective-dispersion  equation  [Equation 
(2.3)],  with  the  retardation  adjusted  to  account  for  mobile  zone  sorption.  Solute  transfer 
between  the  mobile  and  immobile  regions  can  be  accounted  for  in  three  ways  [Rao  and 
others,  1980:1139-1140;  Parker  and  Valocchi,  1986:399;  Brusseau  and  Rao,  1989d:47]: 

(1)  explicitly,  with  a  bicontinuum  approach  describing  the  interchange  between  the 
mobile  and  immobile  zones  governed  by  Ficks  law  of  diffusive  transfer, 

(2)  explicitly,  with  a  bicontinuum  approach  with  interchange  described  by  an 
empirical  first-order  kinetic  expression  (this  assumes  that  the  immobile  region  is  a 
completely  mixed  zone), 

(3)  implicitly,  with  a  lumped  dispersion  coefficient  that  includes  the  effects  of 
sink/source  diffusion  as  a  parameter  which  replaces  the  hydrodynamic  dispersion 
coefficient  (D)  in  the  advection  dispersion  equation  [Equation  (2.3)]. 

First-Order  Rate  Model 

The  first-order  rate  model  describes  a  two-zone  porous  medium  in  which  complete 
mixing  is  assumed  within  the  immobile  region.  This  simplified  model  uses  a  linear 
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expression  as  the  driving  force  approximation.  Solute  transfer  is  assumed  to  be  a  function 
of  solute  concentration  differences  between  the  mobile  and  immobile  zones  and  is  lumped 
into  an  empirical  mass-transfer  term  to  account  for  this  gradient.  These  zones  act  as 
diffusion  sources/sinks  for  solute  flowing  in  the  mobile  regions  [Goltz  and  Roberts, 
1986b:81].  The  model  uses  the  following  equations  [Van  Genuchten  and  Wierenga, 
1976:474]: 
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where  0m  and  0^  are  the  fractions  of  the  soil  filled  with  mobile  and  stagnant  water, 
respectively,  Cm  and  Cjm  are  the  concentrations  in  both  the  mobile  and  stagnant 
regions,  vm  is  the  average  pore-water  velocity  in  the  mobile  liquid,  and  a  is  the  first 
order,  mass  transfer  coefficient. 


This  model  has  been  refined  assuming  linear,  reversible,  equilibrium  sorption  for 
solute  within  immobile  zones,  as  well  as  for  the  "mobile"  solute  and  performing  a  mass 
c  ilance  for  solute  in  the  mobile  region  by  others  [Goltz  and  Roberts,  1 986b:8 1 ;  Brusseau 
and  Rao,  1989d:57;  Harmon  and  others,  1989:412]  to: 
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where  Dm  is  the  dispersion  coefficient  for  the  mobile  zone,  Rm  and  R^  are  the 
retardation  factors  in  the  mobile  and  immobile  zones,  respectively,  defined  as: 
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where  f  is  a  certain  fraction  of  sorption  sites  that  are  in  direct  contact  with  mobile  water. 

Returning  to  the  discussion  concerning  the  first  order  rate  expression  where  we 
make  the  assumption  that  immobile  water  zones  are  perfectly  mixed.  Equation  (2.7) 
becomes  [Goltz  and  Roberts,  1986b:82;  Harmon  and  others,  1989:412;  Brusseau  and 
Rao,  1989d:56;  Brogan,  1991:7]: 

dC 

=  (2.11) 


A  number  of  authors  have  shown  that  the  first-order  rate  model  closely 
approximates  the  tailing  and  asymmetric  affects  expected  of  non-equilibrium  sorption  in 
column  breakthrough  curves  for  certain  parameters  and  conditions  [Van  Genuchten  and 
Wierenga,  1976:479;  Goltz  and  Roberts,  1986b: 85;  Van  Genuchten  and  others,  1977:283; 
Brusseau  and  Rao,  1989d:55-56].  This  model,  like  the  first-order,  kinetic  non-equilibrium 
model  which  follows  below,  avoids  dealing  with  pore  geometry,  which  the  more 
mathematically  rigorous  diffusion  model  uses  with  Ficks  law  of  diffusion.  Due  to  its 
simplicity,  the  first  order  rate  model  has  been  extensively  used  for  modeling  of  physical 
non-equilibrium  processes  in  the  literature. 


2-15 


Diffusion  Model  With  Ficks  Law 

The  second  type  of  physical  non-equilibrium  model  is  known  as  the  diffusion 
model  and  incorporates  Ficks  law  of  diffusion  rather  than  the  empirical  first-order 
formulation  in  Equations  (2.8)  or  (2.1 1)  [Brogan,  1991:8].  Molecular  diffusion  generally 
dominates  microscopic  mass  transfer  processes  and  in  the  liquid  phase  is  controlled  by 
Fickian  (random)  motion  [Weber  and  others,  1991:513].  Several  component  mechanisms 
are  involved  in  the  solute-transfer  process  [Van  Genuchten  and  Wierenga,  1976:473-474; 
Szecsody,  1988:2;  Brusseau,  1989a:6]: 

(1)  advective-dispersive  transport  from  bulk  solution  into  the  boundary  layer  of 
adsorbed  water  surrounding  the  soil  matrix, 

(2)  diffusive  transport  across  the  adsorbed  water  (film  diffusion), 

(3)  pore  and/or  surface  diffusion  within  the  immobile  region  (intra-aggregate 
diffusion  and  within  dead-end  pores). 

The  mobile  phase  is  generally  assumed  to  be  well-mixed,  thus  minimizing  the 
importance  of  intra-aggregate  resistance  within  the  mobile  regions.  All  three  components 
may  be  a  rate-limiting  step.  Intra-aggregate  diffusion  is  usually  the  transfer  process 
designated  as  the  rate-limiting  step  in  sorption  [Brusseau  and  Rao,  1989d:47].  Other 
authors  [Crittenden  and  others,  1986:282]  have  demonstrated  that  intra-aggregate 
diffusion  can  be  modeled  either  as  pore  or  surface  diffusion  and  that  characteristic  skewing 
of  breakthrough  curves  for  tailing  matched  experimental  data. 
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For  diffusion  models,  the  advective-dispersion  equation  is  still  applicable. 
However,  since  this  model  allows  for  concentration  gradients  within  the  immobile  zone 
(vice  use  of  a  perfectly  mixed  zone  assumption  of  the  first-order  model),  the  immobile 
zone  concentration  value  will  represent  a  volume  averaged  solute  concentration  from 
within  the  immobile  zone.  The  immobile  zone  is  generally  represented  by  a  sphere  [Goltz 
and  Roberts,  1986c:  1 140].  The  general  equation  describing  transport  in  the  mobile  zone 
is  [Van  Genuchten,  1985:514;  Goltz,  1986a:8;  Goltz  and  Roberts,  1986b:81-83;  Goltz 
and  Roberts,  1986c:  1 140-1 141;  Goltz  and  Roberts,  1988:41]: 
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and  for  diffusion  into  or  out  of  spherical  aggregates  by  [Parker  and  Vaiocchi,  1986:400]: 


dCo  =  De  d  rr2 
dt  r2R.  3 r  dr 
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giving  the  mean  resident  concentration  in  the  spheres  by: 

Cim=(3/a3)olar2Cadr  (2.14) 

where  0m  and  0;m  are  the  volumetric  water  contents  of  the  inter-aggregate 
(macropore)  and  intra-aggregate  (micropore)  liquid  phases,  Rm  and  R^  are  the 
retardation  factors  of  the  two  regions,  Cm  and  are  the  average  solution 
concentrations  of  the  inter-  and  intra-aggregate  liquid  phases,  respectively,  Dm  is  the 
dispersion  coefficient  for  transport  through  the  macropore  region,  r  is  the  radial 
coordinate  within  the  immobile  zone,  a  is  the  immobile  region  radius  or  layered  half- 
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width,  vm  is  the  average  pore-water  velocity  of  the  macropore  liquid  phase,  z  is  soil 
depth  and  t  is  time. 


The  diffusion  coefficient  for  the  immobile  zone  is  given  by  [Wu  and  Gschwend, 
1988:1373]: 


D  n2 

D,  = - - - 

(l-n)p  +n 


(2.15) 


where  n  is  the  intra-aggregate  porosity,  ps  is  the  dry  solid  density,  D0  is  the  liquid 
diffusion  coefficient,  Dm  is  the  molecular  dispersion  coefficient,  Kp  is  a  partition 
coefficient  and  x  is  tortuosity. 


Within  certain  constraints  the  first-order  model  has  been  shown  to  simulate 
spherical  diffusion  predictions,  but  not  for  all  matrix  geometries  and  non  dibrium 
conditions  [Rao  and  others,  1980: 1 145;  Van  Genuchten,  1985:523;  Parker  and  Valocchi, 
1986:407].  In  fact,  it  may  over-predict  contaminant  mass  remaining  within  the  matrix  at 
high  pumping  rates  [Goltz  and  Oxley,  1991a:554].  Thorbjamarson  [1990:148]  postulated 
that  the  first-order,  mass  transfer  model  may  not  accurately  simulate  the  diffusion  process 
in  situations  of  high  flow  (due  to  pore  water  velocity  influence  on  the  mass  transfer 
coefficient  [a])  or  small  sampling  scale  relative  to  the  diffusion  rate.  Further,  the  first- 
order  model  relies  on  curve  fitting  to  determine  the  mass  transfer  coefficient  at  varying 
velocities  and  is  less  amenable  to  be  used  in  a  predictive  mode  [Goltz  and  Roberts, 
1986b:91].  Finally,  the  spherical  diffusion  model  is  more  firmly  founded  on  a  physical 
basis  and  is  more  amenable  to  independent  determination  of  parameters  [Rao  and  others, 
1980: 1 145;  Goltz  and  Roberts,  1986b:89]. 


2-18 


Pore  Geometry 

As  stated  above,  a  critical  factor  to  the  diffusion  model  is  that  it  requires 
knowledge  of  the  specific  pore  geometry  within  the  immobile  regions  of  the  intra¬ 
aggregate  matrix.  The  equations  above  assume  intra-particle  diffusion  is  occurring  within 
an  ideal  spherical  aggregate.  Several  authors  have  derived  equivalent  equations  for 
diffusion  within  cylindrical  [Van  Genuchten  and  others,  1984:336;  Rasmuson,  1985: 1117], 
rectangular  [Goltz,  1986a:  12-13],  or  layered  shape  geometry’s  [Sudicky  and  Frind, 

1982: 1634;  Rasmuson,  1985: 1 1 17]. 

True  porous  media,  however,  may  be  far  from  this  idealization.  In  an  attempt  to 
apply  the  simplified  model,  Rao  and  others  [1982:348],  demonstrated  that  non-spherical 
aggregates  can  be  represented  by  equivalent  spheres.  Additionally,  a  non-homogenous 
matrix  of  spheres  could  be  characterized  by  a  volumetric,  weighted  average  for  the 
spherical  radius.  These  would  have  a  nominal  radius  which  gives  a  volume  equivalent  to 
that  of  the  aggregate  it  represents  [Crittenden  and  others,  1986:272].  Van  Genuchten 
[1985:520]  extended  the  physical  diffusion  model  to  more  general  conditions  involving 
aggregates  consisting  of  plane  sheets,  solid  cylinders  and  rectangular  prisms.  He 
tabulated  shape  factors  for  these  geometries,  which  could  be  used  with  analytical  solutions 
of  the  advection-dispersion  equation  as  modified  for  spherical  diffusion. 

Using  chemical  engineering  techniques,  another  group  of  authors  proposed  a 
model  to  account  for  differing  aggregate  distributions  which  used  a  swelling  crystal  to 
account  for  size  distribution  within  the  matrix  [Moharir  and  others,  1980: 1795].  The 
effects  of  aggregate  size  distributions  may  be  important  at  the  field  scale,  however,  other 
authors  have  observed  that  its  affects  may  be  ignored  when  the  size  distribution  is 
sufficiently  narrow  [Brusseau  and  Rao,  1989d:55]. 
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Non-Eouilibrium  Model 


Empirical  first-order  models  have  been  discussed  by  several  authors  [Lapidus  and 
Amundson  (1952);  Cameron  and  Klute  (1977);  Rao  and  others  (1985);  Valocchi  (1985); 


Brusseau  and  Rao  (1989b);  Brusseau  and  others  (1991b);  Brogan  (1991);  Haggerty 
(1992)].  In  these  models,  non-equilibrium  is  assumed  to  be  the  result  of  a  time-dependent 
sorption  reaction  at  the  sorbing  surface  [Brusseau  and  Rao,  1989d:43]. 


Cameron  and  Klute  [1977]  published  a  two-site  model  which  assumed  sorption 
sites  were  divided  into  two  fractions,  where  one  has  adsorption  occurring  instantaneously 
and  on  the  othe.  fraction,  it  occurs  in  a  time  dependent  manner.  This  combined  the  first- 
order  chemical  kinetic  model  with  the  linear  equilibrium  model.  The  equations  for  the 
sorbed  component  are  written  as  [Cameron  and  Klute,  1977:183;  Brusseau  and  Rao, 
1989d:43;  Brogan,  1991:12-13]: 


-^  =  m-F)KpC-S„] 


(2.17) 


where  a  is  the  mass-transfer  coefficient  for  chemical  two-site  non  equilibrium  model,  and 
where  F  is  the  fraction  of  immobile  zone  sites  occupied  by  Sj . 


Kinetic  non-equilibrium  formulations  have  been  extended  by  other  authors  who 
developed  a  sorption,  kinetically  influenced  term  model  (SKIT)  which  solves  parallel 
derivations  of  LEA  and  kinetic  non-equilibrium  formulations  for  a  given  solute  transport 
problem  [Bahr  and  Rubin,  1987:439].  Others  have  developed  a  dual  process  (two-site) 
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model  to  simulate  the  BTC  of  PCE,  which  coupled  two,  non-equilibrium  mechanisms,  one 
associated  with  physical  non-equilibrium  and  the  other  that  represented  one  of  several 
sorption-related  processes  (e.g.,  chemical  non-equilibrium,  intra-sorbent  non-equilibrium 
diffusion,  spatial  variations  in  sorption)  [Brusseau  and  Rao,  1989c:241],  Like  the  first- 
order  physical,  non-equilibrium  model,  the  first-order,  kinetic  model  does  not  rely  on  an 
actual  description  of  the  porous  mediums  structure  [Brusseau,  1989a:  16],  Nkedi-Kizza 
and  others  [1984:1 123]  have  shown  that  a  two-site,  first-order  kinetic  non-equilibrium 
model  is  equivalent  to  the  physical  first-order  rate  model  which  uses  mass  transfer  and 
linear  sorption  isotherms  for  simulating  local  equilibrium  conditions.  While  the  two  first 
order  rate  methods  have  been  shown  to  be  mathematically  equivalent,  at  local  equilibrium, 
they  have  different  effects  upon  the  conservative  solute  breakthrough  curves 
[Thorbjamarson,  1990: 16].  This  is  particularly  relevant  in  the  case  of  hydrophobic,  but 
non-sorbing  solutes.  Additionally,  this  model  has  not  been  able  to  satisfactorily  predict 
early  breakthrough  at  higher  velocities  of  breakthrough  curves  for  some  column 
experiments  [Van  Genuchten  and  Wierenga,  1976:475;  Wu  and  Gschwend,  1986:717]. 

Multi-Process  Non-Equilibrium  Model 

Breakthrough  asymmetry  may  not  be  properly  simulated  by  only  one  of  the  non- 
equilibrium  processes.  Recent  research  postulated  that  a  new  model  that  explicitly 
accounts  for  multiple  processes  of  non-equilibrium  operating  in  series  or  parallel  can 
approximate  bicontinuum  models  (physical  non-equilibrium)  for  certain  parameters 
[Brusseau  and  others,  1989b:  1971].  This  model  was  published  with  some  errors,  which 
although  corrected  in  later  publication  [Brusseau  and  others,  1991c:657],  have  vet  to 
correct  mathematical  comparisons  with  data  validations  which  appeared  in  the  original 
[Ball  and  others,  1991:654],  Nevertheless,  an  extension  of  this  model  was  used  to  predict 
sorption,  degradation  and  transport  of  a  contaminant  in  a  packed  column  experiment 
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[Brusseau  and  others,  1992a:  180].  However,  previous  experimental  results  have  not 
shown  that  kinetic  processes  contributed  to  the  rate-limiting  mechanism  for  non-ionic, 
low-polarity  organic  chemicals  [Brusseau  and  Rao,  1991a:  141 ;  Szecsody,  1988:7-8]. 
However,  they  may  be  important  for  other  organic  chemicals  such  as  pesticides  [Brusseau 
and  Rao,  1989e:1691;  Brusseau  and  others,  1991a:134],  Finally,  this  MPNE  model  has 
yet  to  be  tested  adequately  to  support  the  author's  assertions  [Ball  and  others,  1991:653]. 

Pump  and  Treat  Systems 

General 

When  a  site  has  been  well  characterized  for  hydrogeologic  conditions  and 
contaminant  concentrations,  alternatives  for  control  and  remediation  can  be  selected  and 
combined  to  provide  an  overall  strategy  for  cleanup  [Charbeneau  and  others,  1992: 101]  . 
There  are  many  alternatives  to  choose  from  in  handling  ground  water  contamination 
problems,  whether  they  be  containment,  excavation,  pump-and-treat,  vacuum  extraction, 
bio-remediation  or  in-situ  treatment  [US.  EPA,  1992: 130;  Canter  and  others,  1987:467- 
470;  US.  EPA,  1990:1;  Gorelick  and  others,  1993:7],  The  selection  amongst  alternative 
methods  is  beyond  the  scope  of  this  discussion.  Other  authors  [Repa  and  others,  1985:4- 
21  to  4-48;  Canter  and  others,  1987:467-492;  Brubaker  and  Stroo,  1992:167-168] 
provide  general  considerations  for  the  selection  of  technologies  current  at  the  time  of  their 
writing. 


Pump-and-treat  ground  water  extraction  is  the  most  commonly  used  remedial 
technology  for  contaminated  aquifers  [Charbeneau  and  others,  1992:4;  US.  EPA, 
1992:136;  US.  EPA,  1989a:E-l],  A  general  survey  of  1 12  sites  where  pump-and-treat 
was  planned  or  operating  identified  the  following  general  conclusions  [US.  EPA,  1989a.E- 
1  to  E-2]: 
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(1)  the  ground  water  extraction  systems  were  generally  effective  in  containing  the 
contaminant  plumes  and  were  most  effective  at  sites  where  underlying  aquifers  have  high 
inter-granular  hydraulic  conductivity  [Repa  and  others,  1985:5-1], 

(2)  significant  removal  of  contaminant  mass  from  the  subsurface  was  achieved, 

(3)  contaminant  concentrations  usually  decrease  most  rapidly  soon  after  the 
initiation  of  extraction.  After  this  initial  reduction,  the  concentrations  tended  to  level  off 
and  progressed  towards  the  clean-up  goal  at  a  slower  rate  than  expected. 

As  previously  stated,  several  physical  phenomena  have  been  identified  that  tend  to 
interfere  with  aquifer  clean-up  performance  of  pump-and-treat  systems  [US  EPA,  1989a:  1- 
2;  US  EPA,  1990:2]: 

(1)  adsorptive  partitioning  of  contaminants  between  the  ground  water  and  the 
aquifer  materials  (which  this  thesis  is  designed  to  model), 

(2)  aquifer  heterogeneity  (I  assume  homogeneity  for  this  study), 

(3)  low  transmissivity  of  the  aquifer  (which  I  assume  is  high  [greater  than  10'^ 
cm/s]  enough  to  enable  selection  of  pump-and-treat  as  a  viable  option), 

(4)  low  permeability  zones  which  trap  contaminants  within  the  aquifer  (which  are 
neglected  in  this  study  due  to  assumption  of  homogeneity). 
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(5)  immobile  non-aqueous  phase  liquids  (NAPLs)  that  may  contribute  to  a  miscible 
contaminant  plume  by  prolonged  dissolution  by  creation  of  a  separate  phase  at  residual 
saturation  (in  this  instance  pump-and-treat  is  rate-limited  by  how  fast  the  NAPL 
components  can  dissolve,  which  may  be  longer  than  is  reasonable  for  treatment  by  this 
method). 

Ground  water  extraction  technology  is  based  on  two  fundamental  assumptions 
[US.  EPA,  1989a:3-l]: 

(1)  a  well  system  can  produce  ground  water  flow  patterns  that  will  permit  the 
wells  to  withdraw  all  of  the  contaminated  ground  water  from  the  aquifer, 

(2)  contaminants  will  come  out  of  the  aquifer  with  the  water.  For  this  study,  I  will 
assume  that  the  contaminant  of  concern  is  an  organic  (Carbon  Tetra-Chloride  [CTET]) 
which  has  a  high  mobility  as  dissolved  species  in  ground  water. 

Continuous  Pumping 

Continuous  pumping  maintains  an  inward  hydraulic  gradient  and  as  such  pumping 
water  containing  dissolved  contaminants  can  be  addressed  using  standard  well  mechanics 
and  capture  zone  theory  [Charbeneau  and  others,  1992: 103;  US.  EPA,  1990:13]. 
Continuous  pumping  may  be  viable  in  certain  instances  where  the  contaminant  of  interest 
doesn't  sorb  highly.  Brogan  [1991:165-166]  conducted  a  study  using  the  kinetic  approach 
to  adsorption  which  specified  ranges  of  parameters  and  well  combinations  which  were 
used  in  determining  when  to  use  pulsed-pumping  or  continuous  pumping  .  Adams  and 
Virmontes  [1993:5-3]  concluded  that  continuous  pumping  was  viable  in  certain  instances 
as  well,  primarily  when  cleanup  achievement  was  evaluated  in  terms  of  total  remediation 
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time.  Most  pump-and-treat  systems  use  continuous  pumping.  However,  a  few  sites  used 
pulsed-pumping  as  a  technique  to  increase  the  efficiency  of  ground  water  extraction  where 
contaminants  were  subject  to  rate-limited  sorption  [US.  EPA,  1989a:5-l  to  5-2].  The 
tailing  effect  of  the  low  solubility,  sorbed  contaminant  in  these  instances  greatly  affected 
the  economic  and  timely  removal  and  renovation  of  ground  water  due  to  their  tendency  to 
slowly  desorb  into  the  ground  water,  thus  requiring  extended  periods  of  pumping  and 
treating  to  attain  desired  levels  of  restoration  [US.  EPA,  1992:136-137].  Figure  2-2, 
illustrates  tailing  during  pump-and-treat  remediation  showing  that  tailing  entails 
significantly  longer  pumping  times  than  that  predicted  by  simple  equilibrium  models: 
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0  1  2 

Water  Filled  Aquifer  Volumes 


Figure  2.2.  Example  of  Tailing  During  Pump-and-Treat  Remediation 
[Schmelling  and  others,  1992:223;  US.  EPA,  1990:5] 


Another  phenomena  which  occurs  after  the  pumps  are  turned  off  is  known  as 
'rebound.'  In  rebound,  contaminant  desorbs  after  the  induced  gradient  returns  to 
equilibrium.  This  causes  contaminant  concentrations  to  rise  to  a  new  equilibrium  level 
between  that  contained  in  the  immobile  regions  and  the  mobile  regions.  The  result  is  an 
increase  in  contaminant  concentrations  after  the  aquifer  was  supposedly  cleaned-up  by 
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continuous  pumping.  Figure  2.3,  illustrates  the  rebound  aspect  with  respect  to 
contaminant  concentration  levels  over  time: 


Figure  2.3.  Conventional  Pump-and-Treat  With  Rebound 

[Gorelick  and  others,  1993:224] 


Pulsed  Pumping 

An  innovation  in  pump-and-treat  technology  designed  to  increase  the  removal  rate 
for  contaminants  is  pulsed  pumping  [US.  EPA,  1992: 137].  The  principle  behind  pulsed 
pumping  is  that  efficiency  in  contaminant  mass  extraction  will  improve  will  improve  by 
allowing  contaminants  located  in  low  permeability  zones  (such  as  dead-end  pores)  to 
diffuse  outward  into  areas  of  higher  hydraulic  conductivity  during  non-pumping  periods 
[Hall  and  Johnson,  1992:221],  In  this  instance  the  rate  of  mass  removal  is  primarily 
controlled  by  the  release  of  residual  contaminants,  rather  than  by  the  velocity  of  ground 
water  flow.  Thus,  the  pumping  plan  should  accommodate  optimized  "rest"  periods  to 
allow  desorption  from  the  immobile  region  and  terminating  when  "equilibrium"  occurs 
with  the  mobile  region.  Then,  when  extraction  resumes,  the  minimum  volume  of  ground 
water  can  be  removed  at  the  maximum  possible  concentration  for  the  most  efficient 
treatment  [US.  EPA,  1990:25].  Figure  2.4,  illustrates  the  cyclic  nature  of  pulsed  pumping 
and  the  removal  of  rebounded  concentrations  of  residual  contaminants: 
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Figure  2.4.  Pulsed  Pumping  of  Residuals 
[Gorelick  and  others,  1993:225] 

This  increased  efficiency  was  seen  with  a  finite  element  model  [Huso,  1989:4-12], 
and  with  a  kinetic  mathematical  model  [Brogan,  1991:131].  In  the  former  case,  using 
lower  pumping  rate  periods,  increased  average  concentrations  of  organic  species  were 
removed  during  the  "pumping"  periods,  thus  decreasing  the  total  volume  of  water 
requiring  treatment.  In  the  latter  case,  optimization  of  the  SUTRA  code  delineated 
parameters  where  pulsed-pumping  would  be  the  optimal  technology  for  certain  instances 
using  a  kinetic  sorption  model.  Similarly,  using  an  analytical  code  for  a  first-order,  mass- 
transfer,  diffusion  model,  regions  where  pulsed  pumping  would  yield  better  results  than 
continuous  pumping  were  calculated  [Adams  and  Virmontes,  1993:4-19].  In  conclusion, 
an  efficient  pulsed-pumping  alternative  (to  continuous  pumping)  should  reduce  both 
treatment  time  and  volume  pumped  (Powers  and  others,  1991:474). 

Summary 

The  principle  measures  of  operational  effectiveness  of  a  remediation  design  are  the 
general  degree  of  hydrodynamic  control  exerted  and/or  the  general  degree  of  contaminant 
cleanup  achieved.  This  effectiveness  can  be  accomplished  by  examining  the  spatial 
uniformity  of  the  control  exerted  and/or  the  cleanup  level  achieved  [Gorelick  and  others. 
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1993:258-259].  Measures  of  operational  efficiency  for  selecting  remediation  procedures 
are  the  minimization  of  total  costs  required  to  reach  and  maintain  remediation  targets; 
maximization  of  contaminant  removal  per  unit  volume  of  ground  water  pumped  and 
treated;  minimization  of  the  total  volume  of  pumped  and  treated  ground  water;  and  time 
for  completion  of  remediation  to  a  given  water-quality  level  to  which  to  restore  the 
aquifer,  as  well  as  the  most  economical  technology  available  to  reach  that  level. 

Pump-and-treat  systems  are  well  established  as  remedial  techniques  for  restoration 
of  ground  water.  However,  many  systems  have  failed  to  achieve  the  optimal  removal 
efficiencies  due  to  an  incomplete  understanding  of  processes  taking  place  within  the  soil 
matrix  or  overly  optimistic  assumptions  made  in  the  modeling  process.  In  many  cases, 
sorption  played  a  major  part  in  delaying  the  restoration  of  ground  water  quality.  Pump- 
and-treat  technology  is  particularly  good  in  instances  where  the  contaminant  of  concern 
does  not  sorb  strongly  to  the  matrix,  yet  is  hydrophobic.  In  this  instance  high-pumping  at 
a  continuous  rate  would  best  achieve  removal  action,  when  other  alternatives  aren't 
available.  However,  in  the  instance  of  non-equilibrium  sorption,  where  prolonged  periods 
of  pumping  may  occur,  innovations  such  as  pulsed-pumping  can  improve  the  efficiency  of 
contaminant  removal  by  allowing  contaminants  to  reach  equilibrium  during  rest  cycle 
periods  and  consequent  greater  removal  of  contaminant  per  volume  of  water  treated 
[Keely,  1989:9], 

A  number  of  authors  have  voiced  concern  about  the  extensive  use  of  pump-and- 
treat  systems  in  remediation  [Cartwright,  1991:64;  Haley  and  others,  1991:1 19;  Hall  and 
Johnson,  1992:215;  Bredehoeft,  1994:95].  The  basic  problem  is  application  of  this 
technology  to  situations  where  heterogeneities  within  the  sub-surface  interfered  with 
effective  pump-and-treat  capture  of  contaminant  and  sorption  effects  which  extended  the 
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pumping  times  to  uneconomical  proportions  [US.  Congress,  OTA,  1989:151],  Some 
authors  have  suggested  that  containment  may  be  the  optimal  solution  in  certain  situations 
[Mackay  and  Cherry,  1989:630;  Cartwright,  1991:66;  Ross,  1993:92]  or  reliance  on  in- 
situ,  bio-remediation  [Brubaker  and  Stroo,  1992:167-168].  Still  other  authors  have 
sought  to  optimize  the  design  through  use  of  linear  sorption  effects  while  using  monte- 
carlo  or  numerical  models  which  evaluate  uncertainty  in  hydraulic  conductivity  [Ahlfeld, 
1991:848;  Kuo  and  others,  1992:95,105].  Finally,  others  have  used  kinetic  approaches 
towards  optimizing  the  design  [Brogan,  1991;  Haggerty,  1992].  This  thesis  seeks  to  use 
physical  non-equilibrium,  diffusion  oriented  methods  to  analyze  pumping  rates  and 
strategies  to  successfully  remediate  a  hypothetical  aquifer. 

A  further  refinement  to  pump-and-treat,  which  is  beyond  the  scope  of  this  research 
is  combination  of  pump-and-treat  extraction  methods  with  alternative  technologies  such  as 
soil  flushing,  soil  venting,  chemical  enhancements  to  pump-and-treat  or  bio-remediation. 
Discussions  concerning  these  technologies  can  be  found  in  Ross  (1993),  Brubaker  and 
Stroo  (1992),  Bohm  (1992), ,  NGWA  (1992),  Palmer  and  Fish  (1992),  Chambers  and 
others  (1990),  and  Canter  and  others  (1987).  Recent  research  has  focused  on  multi¬ 
process,  non-equilibrium  models  [Brusseau  and  others,  1992a:  175;  Brusseau  (1989b); 
Lewis  and  others,  1987:81;  Jennings  and  others,  1984:120;  Cameron  and  Klute, 

1977:183],  multi-phase  contaminant  transport  (DNAPL's  and  LNAPL's)  [Mayer  and 
Miller,  1990:217]  and  transport  processes  through  fractured  media  [Bear  and  others 
(1993);  Huyakom  and  others,  1983b:841].  The  fractured  media  models  have  great 
potential  for  modeling  the  rapid  movements  through  fractured  soils  or  cracks  in  confining 
beds  which  allow  transmissivity  to  other  aquifers.  These  may  warrant  further  study  and 
development  through  incorporating  physical  non-equilibrium  diffusion  equations  for 
ground  water  transport  problems. 
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I1L  SUTRA  Model  Modification 


Introduction 

The  primary  objective  of  this  section  is  to  develop  a  modification  to  the  USGS 
Saturated-Unsaturated  Transport  (SUTRA)  code  which  describes  contaminant  transport 
by  means  of  advection,  dispersion,  and  physical  non-equilibrium  sorption  (both  first-  and 
second-order)  in  an  aquifer.  Methods  used  will  extend  the  work  of  Goltz  and  Oxley 
[1991a],  Huso  [1989]  and  Huyakom,  et.  al  [1983b].  This  modification  will  be  validated  in 
Chapter  IV  through  comparison  with  an  analytical  code  published  by  Adam's  and 
Viramontes  [1993]  for  a  radial,  pulsed  pumping  simulation.  The  equations,  solutions, 
notation  and  source  code  will  be  based  on  the  expressions  of  sorption  presented  in  their 
respective  papers  as  well  as  in  Bear  and  Verruijt  [1987]. 

This  chapter  begins  with  a  discussion  on  the  characteristics  of  SUTRA.  Then,  the 
governing  equations  of  ground  water  flow  and  contaminant  transport  will  be  presented. 
This  equation  development  will  be  followed  with  a  discussion  of  the  modifications  to 
SUTRA  and  end  with  the  iterative  solution  scheme  used  to  solve  for  contaminant  values 
in  the  modification. 

SUTRA 

General 

SUTRA  is  a  computer  program  which  simulates  fluid  movement  and  the  transport 
of  dissolved  substances  in  a  subsurface  environment.  The  model  employs  a  tw  o- 
dimensional  finite-element  spatial  discretization  and  integrated  finite-difference  time¬ 
stepping  to  approximate  the  governing  equations  that  describe  transport  in  the  ground 


3-1 


water  [Voss,  1984:3].  Transport  of  solute  is  modeled  by  advection,  dispersion  and 
molecular  diffusion  and  assumes  equilibrium  adsorption  on  the  porous  matrix.  It 
provides  a  choice  of  linear,  'Freundlich'  or  'Langmuir'  equilibrium  sorption  models  to 
accomplish  this  equilibrium  assumption.  SUTRA  is  a  public  domain  code  and  one  of  the 
few  which  models  transport  in  both  the  vadose  and  saturated  zones  of  the  sub-surface.  Its 
modular  design  allows  relatively,  straight-forward  modifications  to  account  for  non¬ 
equilibrium  sorption.  This  code  was  primarily  intended  for  two-dimensional  simulation 
of  transport  of  either  solute  or  energy  in  saturated  variable-density  systems  [Voss, 
1984:5].  Models  can  be  extended  into  three  dimensions  by  assuming  a  fixed  parameter 
for  the  aquifer  thickness.  Brogan  has  previously  modified  the  code  to  account  for  first- 
order  kinetic  (chemical)  reactions.  This  thesis  will  modify  the  code  to  account  for  either 
first-  or  second-order  (Fickian)  diffusion  (physical  non-equilibrium)  within  the  matrix 
and  its  influence  on  rate-limited  sorption  behavior  of  a  pump  and  treat  system. 

SUTRA  Numerical  Methods 

SUTRA  simulation  is  based  on  a  hybridization  of  finite-element  and  integrated- 
finite-difference  methods  employed  in  the  framework  of  the  Galerkin  method  of  weighted 
residuals  [Voss,  1984:8],  Standard  finite-element  approximations  are  employed  only  for 
terms  in  the  balance  equations  which  describe  the  fluxes  of  fluid  mass,  solute  mass  and 
energy.  All  other  non-flux  terms  are  approximated  with  a  finite-element  mesh  version  of 
the  integrated-fmite-difference  method.  Time  discretization  in  SUTRA  is  accomplished 
through  implicit,  backward  Euler  finite  differencing. 

Physical-Mathematical  Basis  of  SUTRA  Simulation 

SUTRA  simulation  combines  two  physical  models,  one  to  simulate  the  flow  of 
ground  water  and  the  second  to  simulate  the  movement  of  a  single  solute  in  the  ground 
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water  [Voss,  1984:15].  Haggerty  [1992:15]  modified  the  code  to  account  for  multi¬ 
species  transport  using  the  first-order,  kinetic  algorithms  developed  by  Brogan.  This  was 
accomplished  through  simultaneous  solution  of  mass-balance  equations  at  each  node  and 
time  step  for  the  various  component  reactions  within  the  matrix.  For  the  modification 
studied  in  this  thesis,  a  single-species  solute  will  be  studied.  A  comparison  will  be  made 
for  the  effects  of  first-order  and  second-order  physical  non-equilibrium  models  of 
sorption/desorption  between  the  mobile  and  immobile  regions  assumed  within  the  soil 
matrix  versus  that  predicted  by  an  LEA  model  (SUTRA  unmodified). 

Water  Mass  Balance  Governing  Equations 

The  common  starting  point  in  the  development  of  differential  equations  to 
describe  subsurface  ground  water  movement  and  the  transport  of  solutes  in  porous 
materials  is  to  consider  the  flux  of  water  mass  and  solute  mass  into  and  out  of  a  fixed 
elemental  volume  within  the  flow  domain  [Freeze  and  Cherry,  1979:388].  A 
conservation  of  mass  statement  for  such  an  elemental  volume  is: 
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To  translate  this  statement  into  mathematical  form,  consider  the  volume  flux  into 
and  out  of  a  small  elemental  volume  in  the  porous  medium,  as  illustrated  in  Figure  3.1 
below: 
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Figure  3.1.  Mass  Balance  in  a  Cubic  Element  in  Space 


where  n  is  the  bulk  porosity  of  the  soil  matrix;  Sw  is  the  water  saturation  index  for  the 
soil;  pw  is  the  density  of  water;  and  v  is  the  average  fluid  velocity.  Thus  the  mass  flux, 
using  the  previous  definitions  is  nSw pwv  .  Porosity  may  be  found  as  [Freeze  and 
Cherry,  1979:36;  Voss,  1984:19]: 


where  Vv  =  the  volume  of  the  voids  and  VT  =  the  total  volume  of  the  matrix. 


The  water  saturation  index  (Sw)  for  the  soil  indicates  the  degree  of  saturation  of 
the  soil  matrix.  When  Sw  =  1,  the  void  space  is  completely  filled  with  fluid  and  is  said 
to  be  saturated  [Voss,  1984:19].  When  Sw  <  1,  the  void  space  is  only  partly  water-filled 
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and  is  referred  to  as  being  unsaturated.  In  this  state,  water  adheres  to  the  surface  of  solid 
grains  by  surface  tension  effects  and  the  fluid  pressure  is  less  than  atmospheric.  The 
water  saturation  index  may  be  calculated  by  the  following  equation: 


where  Vw  is  the  volume  of  water  within  the  matrix. 

The  average  fluid  velocity  ( v),  is  the  average  linear  velocity  of  fluid  with  respect 
to  the  stationary  solid  matrix  within  the  soil.  The  average  linear  velocity  does  not 
represent  the  average  velocity  of  the  water  particles  traveling  through  the  pore  spaces,  but 
is  a  macroscopic  term  representative  of  flow  through  the  entire  matrix.  It  may  be  found 
from  the  Darcy’s  Law  equation  as: 

q  K  kkr 

V  =  -  = - VA : =  -pwgVz)  (3.4) 

n  n  nSw\x„ 

where  q  is  the  Darcy  velocity;  K  is  the  hydraulic  conductivity  of  the  matrix;  k  is  the 
solid  matrix  permeability;  kr  is  the  relative  permeability  to  fluid  flow;  pw,  pw,  and  pw 
denote  water's  density,  viscosity  and  pressure  respectively;  z  denotes  elevation;  and 
(Vpw-pwgVz)  =  Vh  is  the  hydraulic  gradient  (head)  for  the  matrix  [Freeze  and  Cherry, 
1979:71;  Voss,  1984:25-27;  Bear  and  Verruijt,  1987:40-41], 

Hydraulic  conductivity  (K)  is  not  a  fundamental  parameter  describing  flow 
through  the  matrix.  It  is  dependent  on  temperature  (which  is  assumed  to  be  constant  in 
the  SUTRA  model)  and  to  a  lesser  extent  on  concentration.  Intrinsic  permeability  (k)  is 
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used  to  describe  flow  in  the  SUTRA  model  [Voss,  1984:37]  as  it  is  independent  of 
pressure,  temperature  and  concentration  (being  a  function  of  the  medium  only  [Freeze 
and  Cherry,  1979:27]).  The  hydraulic  conductivity  may  be  found  from  the  following 
equation: 


where  k  is  the  intrinsic  (specific)  permeability;  p  is  the  dynamic  viscosity  of  the  fluid; 
and  y  is  the  weight  density  (specific  weight)  of  the  fluid  [  Freeze  and  Cherry,  1979:526], 

Using  the  previous  definitions  and  relationships,  we  return  to  the  element  in 
[Figure  3.1]  to  develop  an  equation  for  mass  balance.  SUTRA  conducts  its  flow 
simulation  by  calculating  how  the  amount  of  fluid  mass  contained  within  the  void  spaces 
of  the  solid  matrix  changes  with  time.  The  total  fluid  mass  present  at  any  time  may  be 
given  by: 


Mass  present  =  0pVT 


(3.6) 


where  0  is  the  moisture  content  of  the  soil.  The  volumetric  moisture  content  (0)  of  the 
soil  is  found  by  the  equation  [Freeze  and  Cherry,  1979:39]: 


The  law  of  conservation  of  mass  for  transient  flow  in  a  saturated  porous  medium 
requires  that  the  net  rate  of  fluid  mass  flow  into  any  elemental  control  volume  plus  the 
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net  rate  of  change  of  Ifuid  mass  due  to  reactions  or  sources/sinks  [Figure  3.1]  is  equal  to 
the  time  rate  of  change  of  fluid  mass  storage  within  the  element  [Freeze  and  Cherry, 
1979:64].  In  a  saturated  media,  Sw  =  1  so  the  mass  of  water  in  a  representative  volume  is 
(npw )  and  the  change  in  water  mass  per  unit  volume  of  porous  media  per  unit  time 
[Figure  3.1]  is  given  by  the  equation  of  continuity  for  the  water  [Freeze  and  Cherry, 
1979:531;  Bear,  1972:197;  Bear  and  Verruijt,  1987:57]: 


d(«pw) 

dt 


=  -V(«p  wv)  =  n 


dp* 

dt 


(3.8) 


Now,  since  the  total  mass  of  fluid  in  the  element  consists  of  water  and  dissolved 
solute,  the  flow  simulation  must  keep  track  of  the  total  fluid  mass  and  the  solute  mass 
contained  at  every  point  in  the  simulated  ground  water  system.  The  fluid's  total  density  is 
the  sum  of  that  attributed  to  water  and  the  contribution  from  solute  dissolved  in  the  water. 
For  solute  transport  the  density  is  given  by  [Voss,  1984:18]: 

p  =  p(Qsp  „-Ac-C„)  (3.9) 


where  p  is  the  fluid  density;  C  is  the  contaminant  aqueous  concentration;  and  pw  is 
the  base  density  of  the  fluid  (usually  that  of  pure  water).  Using  the  vectors  from  Figure 
3.1,  we  apply  the  conservation  of  mass  principle  again,  this  time  to  all  of  the  fluid  mass  in 
the  element  [Voss,  1984:33]: 

=  -V(nSw  pv )  +  Q.  +  T  (3.10) 

* 
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where  Qp  is  the  fluid  mass  source  (water  mass  plus  dissolved  solute  mass);  and  T  is  the 
solute  mass  source  (e.g.,  dissolution  of  solid  matrix  or  desorption)  and 


nS 


w 


(3.11) 


Considering  the  left-hand  side  of  Equation  (3.10),  the  actual  amount  of  total  fluid 
mass  contained  in  the  volume  depends  only  on  fluid  pressure  (p)  and  solute  concentration 
(C)  for  isothermal  flow.  The  change  in  total  mass  of  the  fluid  with  time,  assuming  the 
volume  is  constant  is  [Voss,  1984:22;  Bear  and  Verruijt,  1987:57]: 


d(nSw  p)  3 


3(0  p)  dp  3(0  p)  3C 


=  -[(0p)Vr]  =  Vr[ 

3/  dt  dp  dt  dC  dt 


] 


(3.12) 


Saturation  (Sw)  depends  only  on  fluid  pressure  and  porosity  (n)  does  not  depend 
on  concentration.  Since  Sw  is  a  function  only  of  p: 

(313 

dp  dp  dp  dp 


and 


3(0p)  dQiSjp)  nS  5p 
3C  dC  "  3C 


(3.14) 


Substituting  Equations  (3.13)  and  (3.14)  into  Equation  (3.12)  we  obtain  an  expression  for 
the  rate  of  fluid  mass  stored  in  the  representative  elemental  volume: 
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(3.15) 


3Vr(«s.p)1  W  £«£, 

3/  w  dp  dp  dt  dc  dt 

Beginning  with  Equation  (3.8),  I  will  now  discuss  the  effects  of  vertical 
compressibility  on  the  deformable  (consolidating)  porous  medium.  The  specific  mass 
storativity  (S^)  is  the  mass  of  fluid  released  (or  added)  to  a  unit  volume  of  porous 
medium  per  unit  decline  (or  increase)  in  pressure  [Bear,  1972:204].  It  is  related  to 
pressure  changes  by: 

^Sw  =  -p[p«  +  a'(l-»)]^  =  pVv  (3.16) 

where  p  is  the  coefficient  of  compressibility  of  water;  and  a'  is  the  soil  coefficient  of 
compressibility.  Aquifer  storativity  under  fully  saturated  conditions  is  related  to  the 
factor,  d(np)/dp  as  follows  [Voss,  1984:23]: 

^  =  P (3.17) 

dp 


Expanding  Equation  (3.17)  we  obtain  the  relationship  [Bear,  1972:205;  Voss,  1984:23]: 


d(«p) 

dp 


=  ps0p  =  p 


dn  dp 
dp  dp 


(3.18) 


Substituting  Equation  (3.18)  into  Equation  (3.15)  and  then  substituting  the  resultant 
equation  into  Equation  (3.10),  in  addition  to  using  Equation  (3.4)  for  the  Darcy  velocity 
we  obtain  the  more  exact  form  of  the  fluid  mass  equation  which  is  used  in  SUTRA 
(ignoring  T)  [Voss,  1984:34;  Bear  and  Verruijt,  1987:63]: 
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(3.19) 


kk  p  QS  dp  dp 

QP  =  V[(^-)(Vp  -  pgVz)]  +  (Sw  pS0p  +  «P  -~)  ^  +  («S- 

Contaminant  Transport  Equation 

SUTRA  models  solute  mass  transport  through  use  of  an  advection-diffusion 
equation,  with  the  addition  of  a  dispersion  term  to  approximate  the  effects  of  mixing  of 
waters  with  different  concentrations  moving  both  faster  and  slower  than  the  average 
velocity  (v  )[Voss,  1984:38].  This  equation  represents  the  total  mass  of  a  single  species 
of  adsorbent  by  summing  separate  equations  for  the  mass  balance  of  that  stored  in 
solution  and  that  stored  on  the  solid  grains  (adsorbate).  The  derivation  of  the  equations 
for  contaminant  mass  balance  follows  the  same  pattern  as  that  done  above  for  the  ground 
water  which  culminated  in  Equation  (3.19).  The  separate  balances  for  a  single  species 
stored  in  solution  (solute  and  on  the  solid  grains  (adsorbate),  are  expressed,  respectively, 
as  follows  [Voss,  1984:40]: 

=  -/_  V(0pv  O + v[0p(D + d;>vcj + oPr„ + qpc  (3.20) 

dt 

and 


d[(\-n)psF] 

dt 


/  +  (l-n)p,r. 


(3.21) 


where  F  is  the  mass  of  solute  adsorbed  on  the  solid  per  unit  mass  of  solid. 

These  in  turn  are  coupled  additively  to  yield  the  general  form  of  the  SUTRA  integrated 
equation  used  for  total  species  mass  balance  [  Bear  and  Verruijt,  1987:170;  Voss, 
1984:42]: 
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= -  vis.pc^)]  -  + e,p,r, + ePr + q„c  (3.22) 

a/  a/ 

where  0S  is  the  solids  volumetric  fraction,  ps  is  the  solids  density,  Ts  is  the  rate  of 
production  of  the  solute  per  unit  mass  of  solid,  and  QpC*  is  the  dissolved  species  mass 
added  by  a  fluid  source  with  concentration  C*. 


The  last  term  in  Equation  (3.22)  is  represented  by  some  as  two  separate  terms,  and  can  be 
found  by  the  equation  [Bear  and  Verruijt,  1987:170]: 

QpC=-PC+RqCr  (3.23) 


where  P  is  the  rate  of  water  withdrawn  at  a  pumping  well,  Rq  is  the  rate  of  water  added 
at  an  injection  well,  and  CR  is  the  solute  concentration  in  water  injected  into  the 
medium.  The  term  0S  may  be  found  by  the  equation: 


(3-24) 


Advection,  Dispersion  and  Diffusion 

There  are  four  components  taken  into  account  in  the  construction  of  the  solute 
mass  balance  equation  [Equation  (3.22)]  [Bear  and  Verruijt,  1987:168].  The  first  is  the 
quantity  of  the  pollutant  entering  and  leaving  a  control  volume  around  a  considered  point 
by  advection,  dispersion  and  diffusion.  The  total  flux  portion  of  Equation  (3.22)  is 
represented  by  : 

qc.lnlal  =  Cq-qDVC-qDyC  (3.25) 
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where  D  is  the  coefficient  of  mechanical  dispersion;  D*  is  the  coefficient  of  molecular 
diffusion  in  a  porous  medium. 

Advection  is  caused  by  motion  of  the  fluid  into  the  control  volume  [Freeze  and 
Cherry,  1987:389].  This  is  represented  by  the  first  term  in  Equation  (3.25),  where  the  net 
rate  of  solute  inflow  is  -V(0pv).  Mechanical  dispersion  is  mixing  that  occurs  as  a 
consequence  of  local  variations  in  velocity  around  some  mean  velocity  flow  [Domenico 
and  Schwartz,  1990:369].  The  coefficient  of  mechanical  dispersion  (D)  is  the  result  of 
the  motion  of  the  water  in  the  pore  spaces,  due  to  differing  velocities  which  the  water  will 
travel  in  the  pore  space  (channels),  differences  in  surface  area  and  roughness  relative  to 
the  volume  of  water  in  individual  pore  channels  and  tortuosity  related  to  the  branching 
and  inter-fingering  of  pore  channels  [Freeze  and  Cherry,  1979:75-76].  A  number  of 
authors  have  discussed  methods  of  calculating  the  coefficient  of  dispersion  (D)  [Freeze 
and  Cherry,  1979:389-390;  Bear  and  Verruijt,  1987:161-164].  SUTRA  [Voss,  1984:48- 
50]  uses  a  2  x  2  matrix  to  compute  mechanical  dispersion  in  the  manner  suggested  by 
Bear  and  Verruijt.  The  flux  associated  with  molecular  diffusion  ( D'd )  at  the  microscopic 
level  originates  because  of  diffusion  related  to  mixing  caused  by  random  molecular 
motions  caused  by  the  thermal  kinetic  energy  of  the  solute  [Domenico  and  Schwartz, 
1990:367].  Its  effects  are  represented  by  the  third  term  in  Equation  (3.25).  Values  for 
the  molecular  diffusion  coefficient  can  be  calculated  by  the  empirical  methods  of  Bear 
and  Verruijt  [1987:164-165]  or  by  some  more  general  empirical  approaches  to 
approximate  an  effective  diffusion  coefficient  (D'd  )  in  place  of  the  molecular  diffusion 
coefficient  (D‘d)  along  with  some  tabular  values  for  ionic  diffusion  coefficients. 
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Equilibrium  Sorption  Isotherms 

The  next  term  of  Equation  (3.22)  is  that  representing  pollutant  leaving  the  fluid 
phase  through  adsorption  through  the  water-solid  interface  as  a  result  of  chemical  or 
electrical  reactions.  This  is  represented  by  the  volumetric  adsorbate  term  (f)  in  Equations 
(3.20)  and  (3.21),  which  denote  the  quantity  of  solute  that  leaves  the  water  by  such 
mechanisms,  per  unit  volume  of  porous  medium,  per  unit  time  [Bear  and  Verruijt, 
1984:168].  This  represents  adsorption.  Within  SUTRA,  the  volumetric  adsorbate  source 
term  (f)  is  represented  by  equilibrium  isotherms,  which  relate  the  rate  of  transfer  between 
the  mobile  and  immobile  zones  per  those  found  in  Table  2.1.  In  general  these  isotherms 
take  the  form: 

—  =  (Const.)—  (3.26) 

dt  dt 

which  reduces  the  number  of  unknowns  in  Equation  (3.22)  from  two  unknowns  (F  and  C) 
to  one  unknown  (C).  The  three  isotherms  used  by  SUTRA  are  either:  linear  sorption, 
Freundlich  equilibrium  sorption,  or  Langmuir  equilibrium  sorption.  Should  adsorption 
not  be  a  factor,  the  code  makes  provision  to  zero  this  term  out.  Later,  we  will  discuss  the 
methods  used  to  represent  this  term  with  a  second-order  diffusion  equation. 

Sources  and  Sinks  Within  the  Liquid  Phase 

Sources  and  sinks  of  the  solute  are  expressed  by  the  third  and  fourth  terms  in 
Equation  (3.22),  +0pT  .  These  result  from  various  processes,  e.g.,  chemical 
reactions  among  components  within  the  liquid,  radioactive  decay  and  bio-degradation 
[Bear  and  Verruijt,  1987:173;  Voss,  1984:46].  Bear  and  Verruijt  [1987:173]  provide  a 
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generic  discussion  of  sink  terms  and  decay  constants,  as  do  Domenico  and  Schwartz 
[1990:475-477]. 

Solute  Mass  Additions 

The  final  term  of  Equation  (3.22)  is  that  attributed  to  solute  which  may  be  added 
by  injecting  polluted  water  or  through  removal  of  mass  through  pumping  (  QpC*  ). 

This  term  accounts  for  additions  in  mass  due  to  recharge  of  the  aquifer  by  contaminated 
water  by  a  given  concentration,  or  removal,  through  pumping  of  contaminated  water  from 
the  medium.  This  addition,  serves  to  balance  the  equations  used  in  tracking  the  mass 
balance  throughout  the  system. 

Implementing  Non-Equilibrium  Sorption  in  SUTRA 

General 

As  discussed  in  Chapter  II,  physical  non-equilibrium  models  are  commonly  based 
on  a  two-zone  description  of  a  porous  medium,  with  transport  in  the  mobile  zone 
governed  by  the  advective-dispersive  equation  (Equation  [2.3])  with  the  retardation  factor 
adjusted,  as  in  Equation  (2.9),  to  account  for  a  mobile  zone  sorption  [Harmon  and  others, 

1 989:409-4 1 0].  The  first-order  rate  model  (Equation  [2.11])  assumes  that  the  immobile 
zone  is  completely  mixed  and  that  mass  transfer  across  the  mobile-immobile  interface  is 
described  by  a  first-order  expression.  A  second-order  rate  model  (Equations  [2.13]  and 
[2. 14])  is  one  which  considers  the  transfer  of  solute  within  the  immobile  region  as  being 
diffusion  driven,  using  Fick's  law  of  diffusion.  SUTRA's  algorithms  are  structured  for 
local  equilibrium  conditions.  However,  it  does  accommodate  non-equilibrium  sorption 
models,  through  use  of  kinetic  sorption  coefficients  to  formulate  a  value  for  a  specific 
sorption  rate  (f„),  which  in  turn  can  be  used  in  calculation  for  a  value  of  the  volumetric 
adsorbate  source,  (f)  for  use  in  Equations  (3.20)  and  (3.21).  Traditional  approaches 
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couple  Equations  (3.20)  and  (3.21).  Equations  (2.13),  (2.14)  and  (3.22)  are  then  solved 
simultaneously  for  the  spherical  diffusion  model  to  account  for  three-dimensional 
transport  in  the  soil  matrix.  Bibby  [1981 :1075-1076]  demonstrated  use  of  the  three- 
dimensional  equations  above,  while  analytically  solving  for  one-dimensional  diffusion 
within  a  spherical  block  for  Equation  (2.17).  This  implementation  requires  storage  of 
previous  immobile  zone  solutions.  In  contrast,  Huyakom,  et.  al  [1983b:842]  used  a 
finite-element  solution  of  the  immobile  zone  equation  [Equation  (2.16)]  in  a  fractured 
media  model.  The  solution  flexibility  and  smaller  storage  requirements  of  this  method 
favor  the  latter  approach.  Therefore,  transport  in  the  spherical  matrix  block,  which 
represents  immobile  zones,  will  be  approximated  by  a  one-dimensional,  finite-element 
solution  that  is  sequentially  coupled  with  the  two-dimensional  approximation  of  the 
mobile  zone  via  leakage  flux  terms. 

The  conceptual  matrix  models  used  in  this  thesis  effort  are  shown  graphically  in 
Figures  3.2  and  3.3,  below: 


Figure  3.2.  Dual  Porosity  Media-Spherical  Blocks 
[Huyakom  and  others,  1983b:843] 
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Figure  3.3.  Discretization  of  Spherical  Immobile  Region 
[Huyakom  and  others,  1983b:843] 

Huyakom  [1983b:843]  solved  for  the  concentration  values  within  the  spherical 
immobile  region  represented  by  Figure  3.3  by  the  equation: 


dC  D  d  ,  dC 
dr  r2  Sr  dr  ’ 


(3.27) 


and  the  concentration  flux  at  the  boundary  of  the  layers  of  the  immobile  and  mobile  zones 
of  Figure  3.2  by: 


f 


3De  dCim 


a  dr 


(3.28) 


This  expression  simply  relates  the  mass  flux  of  contaminant  flowing  into  or  out  of  the 
immobile  zone  to  the  spatial  gradient  at  the  zone’s  perimeter. 


The  overall  solution  procedure  for  the  non-equilibrium  equations,  either  Equation 
(2.11)  for  first-order  or  Equation  (2.13)  for  second  order,  follows  an  iterative  solution 
scheme  set  forth  by  Huyakom,  et  al  [1983b:846]  for  fractured  media  flow  and  Miller  and 
Ra'oideau  [1993:2229]  for  the  operator-splitting  technique.  Using  the  split-operator 
approach,  the  governing  equations  are  separated  into  transport  (for  the  mobile  zone)  and 
reaction  operators  (diffusion  equations  for  the  immobile  zone)  and  solved  sequentially. 
This  splitting  allows  the  separation  of  the  short-time  scale  process  (mobile  zone)  from  the 
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long-time  scale  process  (immobile).  Operator  splitting  leads  to  smaller  systems  of 
equations,  which  can  be  solved  faster  than  coupling  the  diffusion  equation  into  Equation 
(3.22).  As  outlined  above,  my  approach  will  entailed  solving  Equation  (2.10)  for  first 
order  or  Equations  (3.27)  and  (3.28)  for  second-order  diffusion,  and  finding  the  fluxes 
using  the  split  operator  approach.  The  results  will  be  then  be  sequentially  used  to  solve 
Equation  (3.20)  at  each  time  step,  for  the  mobile  region  concentration  at  the  respective 
nodes.  Miller  and  Rabideau  [1993:2229]  further  discussed  the  numerical  stability  and 
computational  advantages  of  this  method  in  greater  detail. 


First  Order.  Non-Equilibrium  Assumption 

Assuming  that  Ts  =  0  simplifies  Equation  (3.21)  and  the  remaining  term,  f  in 
the  equation  is  the  rate  of  solute  transfer  from  the  matrix  block  to  the  mobile  zone  per 
unit  volume  of  the  matrix  block.  Van  Genuchten  and  Wierenga  [1976:474]  established  a 
first-order  relationship  for  the  immobile  zone  concentration,  which  is  given  by  Equation 
(2.7).  The  concentration  within  the  immobile  region  is  assumed  to  be  a  volumetric 
averaged  constant  and  is  equal  to  the  concentration  at  the  boundary  with  the  mobile  zone. 
The  first  order  rate  constant,  for  layered,  spherical  aggregates,  corresponds  to  the 
following  equation  [Van  Genuchten,  1985:522;  Goltz,  1986a:56]: 


a  = 


22.68 DQ,m 

_ €  Im 


(3.29) 


where  De  is  the  effective  diffusion  coefficient,  9,m  is  the  immobile  zone  porosity  and  a 
is  the  layered  half  width  for  the  immobile  region  (Equation  [3.27]  is  the  same  if  the 
spherical  immobile  region  radius  [a]  is  known  as  well).  This  relationship  is  used,  in  turn, 
to  solve  Equation  (3.21),  where  the  resulting  equation: 
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d[(l-n)p,Fj 

dt 


=  f  =  a(Cm  -Cjm) 


(3.30) 


is  used  to  solve  the  SUTRA  integrated  equation  for  total  species  mass  balance 
(Equation[3.22J). 

Second-  Order.  Non-Equilibrium  Assumption 

Instead  of  using  the  LEA  to  further  simplify  the  solution  to  equations  (3.20)  and 
(3.21),  an  expression  for  f  in  Equation  (3.21)  may  be  found  by  solving  the  one¬ 
dimensional  diffusion  equation  (Equation  [2.13])  for  Ca  using  the  following  boundary 
conditions,  from  Figure  3.3: 

Ca(a)  =  Cm  (3.31) 


and 


dc 


dt 


=  0 


r=0 


(3.32) 


The  non-equilibrium  transport  problem,  Equation  (3.20),  is  solved  through  iterative 
solution  of  Equations  (3.27)  and  (2.28)  using  the  split  operator  technique,  to  determine 
the  mass  flux  from  the  immobile  zone  into  the  mobile  zone,  using  the  boundary 
conditions  of  Equations  (3.31)  and  (3.32). 
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Iterative  Solution  Scheme 

For  second-order  calculations,  at  the  initial  time  step,  the  mobile  and  immobile 
concentrations  are  assumed  to  be  equal  (at  the  boundary  of  the  spheres).  Thus  f  =  0  at 
all  points.  Then,  using  a  modified  form  of  SUTRA,  Equation  (3.21)  is  solved  for  the 
concentration  of  the  mobile  zone  at  Cmk+1  (superscript  represents  time)  iteration  at  each 
node.  The  Cmk+1  values  (where  Cm  is  the  mobile  zone  concentration)  at  each  node  then 
serve  as  a  boundary  condition  [Equation  (3.3 1)]  for  the  immobile  zone  concentration 
(Qmk)  calculation  at  each  node  using  Equation  (3.28).  The  mobile/immobile  zone 
interface  iteration  can  be  summarized  by  the  following  equation: 

c.( a)  =  Cj1"  (3.32) 

Then  the  immobile  zone  concentrations  within  the  spherical-immobile  regions  are 
calculated  using  a  one-dimensional  finite  element  code  [Pepper  and  Heinrich,  1992:226] 
with  Cak  as  an  initial  condition  estimate  for  the  calculation  of  a  revised  value  estimator 
for  Cak+I.  Finally,  the  flux  into  or  out  of  the  spherical  immobile  region  (f)  is  computed 
using  Equation  (3.28)  and  used  as  a  source  term  ( f )  in  the  mobile  zone  calculation  of 
Equation  (3.20).  For  the  next  cycle,  the  improved  estimate  for  the  Cmk+1  values  at  each 
node  are  used  to  update  the  concentration  values  for  the  immobile  zones  in  step  two 
above,  and  conduct  another  series  of  one-dimensional  Galerkin  solutions  in  the  same 
manner  as  before. 


HL  Analysis  and  Evaluation  of  Results 


Introduction 

This  chapter  presents  results  using  modifications  made  to  the  SUTRA  code  based 
on  the  analysis  developed  in  Chapter  HI.  First,  model  verification  is  asserted  by  validating 
the  modified  SUTRA  code  with  an  analytical  code  developed  by  Adams  and  Virmontes 
(1993),  in  a  one-dimensional  (1-D),  radial  pumping  test,  for  both  first-  and  second-order 
diffusion  effects  that  emulate  rate-limited  sorption  and  desorption .  Then,  a  hypothetical 
two-dimensional  (2-D),  contaminated,  sandy  aquifer  was  developed  using  realistic 
parameters  to  facilitate  the  remainder  of  the  analysis.  The  discussion  of  aquifer 
parameters  included  analysis  of  input  parameters  for  placement  within  the  non-equilibrium 
regime.  Repeated  simulations  on  this  hypothetical  aquifer  provided  a  means  of  finding 
discrepancies  between  the  first-  and  second-order,  rate-limited  model  predictions  of 
mobile  zone  concentrations  and  cleanup  times.  These  simulations  also  provided  insight 
into  the  effectiveness  of  pulsed-pumping  as  a  remediation  strategy  when  a  pump  and  treat 
system  is  considered. 


Validation  of  Model  Modifications 

Introduction 

The  basis  for  model  verification  was  to  compare  this  model  with  an  existing  model 
formulated  using  the  same  sorption/desorption  assumptions  and  through  incorporating 
changes  which  had  been  separately  validated  by  their  originators.  In  Chapter  II  a  review 
of  the  literature  related  to  modeling  sorbing  solute  transport  was  conducted.  The  review 
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focused  on  models  incorporating  either  equilibrium  or  rate-limited  sorption,  with  the  rate 
limitation  described  by  either  a  first-order  equation,  or  by  Fickian  diffusion  in  immobile 
regions  of  an  aquifer. 

The  literature  review  discussed  models  presented  by  various  researchers  which  met 
the  above  criteria.  For  purpose  of  validation  comparison,  the  analytical  model  presented 
by  Adams  and  Viimontes  [1993]  was  used.  Their  model  compared  favorably  with  another 
analytical  model  [Goltz  and  Oxley,  1991a]  and  a  numerical  model  [Huso,  1989].  The 
Adams  and  Virmontes  model  solved  a  similar  sorption/desorption  equation  set  to  that  used 
in  SUTRA,  for  a  radial  pumping  simulation,  for  both  first-  and  second-  order  cases.  The 
situation  they  simulated  assumed  a  non-sorbing  solute  (retardation  coefficient[R]  of  1.0)  in 
a  radial  flow  field.  The  simplified  geometry  allowed  the  use  of  1-D  equations  which 
yielded  an  analytical  solution.  Another  unique  feature  of  the  Adams  and  Virmontes 
[1993:4-25]  code  was  that  it  allowed  specification  of  arbitrary  initial  conditions  and 
pulsed-pumping  within  the  aquifer. 

The  separate  portions  of  code  combined  for  the  modified  SUTRA  have  be,  n 
validated  individually.  One  reason  for  the  choice  of  SUTRA  was  its  history  of  validation 
[Voss,  1984: 180];  peer  review  [Van  der  Heijde  and  Elnawawy;  1993:C.4-2-l];  and 
application  to  field  problems  by  many  users  [Kolodny,  1989:147;  Van  der  Heijde  and 
Elnawawy,  1993:C.4-2-l].  The  primary  modification  to  the  SUTRA  code  as  outlined  in 
Chapter  III,  was  modeling  of  the  immobile  region  through  use  of  a  1-D  finite  element 
code.  The  immobile  region  code  was  based  on  a  method  developed  by  Pepper  and 
Heinrich  [  1992:226]  and  the  method  had  been  validated  separately  against  analytical 
solutions  to  diffusion  problems.  Huyakom  et  al  [  1983b:847]  validated  the  approach  used 
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in  the  transport  portion  of  the  code  through  use  of  a  dual-porosity  model  in  solution  of  2- 
D  problems  involving  transport  in  fracture  and  radial  spherical  matrix  diffusion.  Finally, 
split  operator  methods  used  to  solve  the  transport  and  reaction  equations  separately  at 
each  time  step  have  been  shown  to  converge  and  numerical  error  minimized  when  the  size 
of  the  time  step  chosen  was  small  compared  to  the  time  scale  of  the  transport  or  reaction 
process  under  consideration  [Miller  and  Rabideau,  1993:2233-2234].  The  SUTRA 
modification  used  five  (5)  time  steps  in  the  spherical  diffusion  regime  for  every  time  step 
taken  within  the  transport  or  reaction  iterations.  In  all  cases,  the  transport,  reactive  and 
diffusive  module  components  maintained  Courant  numbers  less  than  one 
(  Cr  =  vxAx/At  <  1  )  and  mesh  Peclet  numbers  less  than  two  (  Pe  =  vxAx/D„  <  2  ) 
throughout  the  radial  simulation,  where  Ax  is  the  macroscopic  grid  spacing  and  At  is  the 
time  step.  This  was  particularly  important  as  numerical  instabilities  are  minimized  and 
mole  cular  diffusion  effects  are  ascendant  when  Pe  <  2  and  Cr  <  1  [Goltz  and  Roberts, 
1986b:89;  Anderson  and  Woessner,  1992:327]. 

Radial  l  est  Input  Parameters 

In  this  simulation,  model  comparisons  were  made  between  pulsed  pumped  and 
continuous  pumped  aquifers  for  both  first-order  equations  and  second-order  diffusion 
within  the  immobile  zones.  The  pulsed-pump  simulation  for  the  first-order  simulation 
consisted  of  cycling  the  pump  on  and  off  at  60-day  intervals  for  a  total  duration  of  120- 
days.  There  was  no  initial  hydraulic  gradient  in  the  aquifer,  hence  heads  were  everywhere 
equal.  The  pumping  rate  was  constant  throughout  the  pumping  period.  The  initial  mobile 
and  immobile  contaminant  concentration  distributions  were  set  at  1 .0  (dimensionless) 
throughout  the  56-meter  diameter  site  of  contamination  and  0.0  for  both  mobile  and 
immobile  regions  outside  the  site  boundaries.  A  radial  mesh  was  generated  by  SUTRA 


using  a  plume  radius  (rp)  of  28-meters,  with  mesh  points  generated  at  0.07  meter  intervals 
near  the  pump,  increasing  to  0.7  meter  intervals  at  the  outer  portion  of  the  domain. 
Appendix  A  lists  the  other  input  parameters  used  for  this  simulation. 

Results  of  Validation 

The  first-order  diffusion  modification  to  SUTRA  showed  good  agreement 
throughout  the  simulation.  Similarly,  the  second-order  diffusion  modification  compared 
favorably  at  the  initial  time  steps  for  concentration,  as  shown  in  Table  4.1  below. 
However,  the  second-order  model  showed  significantly  more  extraction  taking  place 
(and  consequently  less  retardation)  at  larger  simulation  times  than  the  analytical  model. 
This  difference  is  attributed  to  a  number  of  factors: 

a.  The  difficulty  in  simulating  a  1-D  flow  field  uniformly  across  the  SUTRA  finite 
element  grid  to  match  the  analytical  1-D  solution. 

b.  Instability  in  the  Laplace  transforms  and  Green’s  functions  of  the  analytical 
model.  As  a  result  of  these  instabilities,  the  analytical  code  was  unable  to  compute  a 
solution  during  the  pump-off  cycle  for  the  second-order  model. 

c.  SUTRA's  boundary  condition  routines  did  not  allow  for  a  complete 
representation  of  the  conditions  used  at  the  outer  boundary  for  the  analytical  code. 

Despite  the  divergence  of  the  second-order  model  solution  at  larger  values  for 
time,  its  solution  was  qualitatively  in  agreement  with  the  analytical  model.  The  SUTRA 
modification  was  considered  valid  for  employment  in  modeling  contaminant  transport  in 
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the  soil  matrix  due  to  its  similar  form,  behavior  and  comparable  results  to  the  first-order 
solution  for  pulsed  pumping. 


TABLE  4.1 

COMPARISON  OF  FIRST-  AND  SECOND-ORDER  NODES  TO  ANALYTICAL 

MODEL 

(VALUES  ARE  CONCENTRATION  AT  THE  PUMPING  WELL) 


lst-Order 

lst-Order 

2nd-Order 

2nd-Order 

Day 

Analytical 

SUTRA 

Analytical 

SUTRA 

Code 

Modified 

Code 

Modified 

24 

0.999 

0.999 

0.999 

0.997 

60 

0.621 

0.655 

0.729 

0.524 

120 

0.620 

0.673 

N/A 

N/A 

Motel  Aquifer  Parameters 
General 

A  search  of  the  literature  was  conducted  to  determine  parameters  which  were 
representative  of  “reality”  (in  the  absence  of  field  data)  for  input  into  the  model  simulation. 
Orientation  was  on  the  well  known  Borden  experiment  [Roberts  and  others,  1986; 

Mackay  and  others,  1986].  This  case  was  chosen  because  it  provided  controlled  data,  had 
been  studied  at  length  and  published  and  cross-referenced  parameters  were  available. 
Additionally,  this  case  would  enable  achievement  of  all  of  the  assumptions  made  in 
Chapter  I.  For  ease  of  solution,  a  homogeneous  aquifer  simplification  was  made,  even 
though  it  has  been  postulated  that  Borden's  retardation  behavior  was  due  more  to 
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heterogeneity  in  the  matrix  than  rate-limited  sorption  [Burr  and  others,  1994:813].  Pulsed 
and  continuous  pumping  at  a  constant  rate  when  the  pump  was  on,  were  employed  at  the 
injection  point  to  simulate  remediation  of  the  mobile  water  within  the  aquifer  to  a 
hypothetical  clean-up  standard.  The  plume  radius  (rp)  was  defined  as  the  radius  at  which 
the  concentration  was  1%  of  the  cleanup  standard  (5  x  10' ^  g/1).  Appendix  B, 
summarizes  geophysical  factors  which  were  used  in  the  simulations. 

Validation  of  Non-Equilibrium  Regime 

Many  authors  have  postulated  regimes  where  non-equilibrium  was  considered 
ascendant  in  so  far  as  its  effects  on  dispersion.  Brusseau  and  Rao  [1989d:43]  stated  that 
in  regions  where  the  Peclet  number  was  less  than  20  (  Pe  <  20  ),  dispersion  would  have 
significant,  if  not  dominant  effects;  whereas  non-equilibrium  effects  would  be  dominant  at 
higher  Peclet  number  regimes.  This  translates  to  non-equilibrium  effects  near  the  pumping 
well  and  dispersive  effects  further  away  from  the  pumping  well.  Computation  of  Peclet 
numbers  for  various  radii  from  the  injection  point  verified  that  Pe  numbers  were  less  than 
20  for  all  points  within  the  plume  radius  from  the  proposed  extraction  point  and  lower  for 
outer  points  (  Pe  <  10  )  within  the  mesh  for  all  immobile  zone  radii  (a),  pumping  rates 
(Qw),  and  diffusion  coefficients  (De)  in  Appendix  B.  Simulation  in  these  regimes 
attempted  to  capture  non-equilibrium  effects  at  the  lead  portion  of  their  applicability. 
Unfortunately,  this  low  Peclet  number  regime  made  the  hydrodynamic  dispersion  effects 
difficult  to  distinguish  from  non-equilibrium  effects. 

Finite  Element  Mesh 

As  outlined  above,  some  authors  had  specified  mesh  selection  using  the  equations 
below  to  eliminate  overshooting  and  reduce  inaccuracy  [Goltz  and  Roberts,  1986b:89j: 
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(4.1) 


D. 


<2 


Ax 


<,2 


(4.2) 


Anderson  and  Woessner  [1992:327]  further  refine  these  criteria  by  stating  that  the  grid 
should  be  designed  so  that: 


A/<4at 

(4.3) 

a  A/ 

A t<  — 

(4.4) 

v 


The  hypothetical  aquifer's  finite-element  mesh  was  designed  by  using  Equations 
(4.1)  through  (4.4),  to  specify  the  maximum  grid  element  size  (Al)  and  time  increment  (At) 
which  would  maintain  the  Pe  <  2  in  the  macropores  for  the  plume.  Therefore,  a  uniform 
finite  element  mesh  28-meters  in  length  and  15-meters  in  width  was  generated  by  SUTRA 
for  model  formulation.  A  maximum  At  value  of  1.0  day  was  used,  although  smaller 
intervals  (0.1  -  0.5  day)  were  preferable  and  used  during  pumping  periods  for  second- 
order  simulations  on  selected  intervals  and  0.01  day  for  first-order  simulations  throughout 
to  improve  accuracy  and  stability.  The  contaminant  plume  was  injected  at  a  point  four 
meters  in  from  the  left  boundary  and  centered  laterally.  The  primary  intention  for 
selecting  the  size  of  the  mesh  was  the  desire  to  reduce  computational  time  required  for 
solution  of  the  second-order  modification.  The  ratio  of  longitudinal  to  transverse 
dispersivity  (12:1)  resulted  in  a  distinctly  elongated  plume  in  the  direction  of  ground  water 
flow,  which  is  a  characteristic  feature  consistent  with  field  observations  often  observed  at 
ground  water  contamination  sites  [Palmer  and  Johnson,  1989:9].  Figure  4.1  illustrates  the 
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rectangular  mesh  generated  and  its  associated  boundary  conditions  for  the  contaminant 
plume: 


Regional  Ground  Water  Flow  Direction 


No  Flux  Boundary 


Constant  Head 
Boundary 


Constant  Head 
Boundary 


No  Flux  Boundary 

Figure  4.1.  Finite  Element  Mesh  and  2-D  Boundary  Conditions  for  Hypothetical  Model 


Immobile  Region  Mesh 

Chapter  III  discussed  modeling  of  the  immobile  region  using  a  finite  element  code 
developed  by  Pepper  ai  d  Hienrich  [1992].  The  finite  element  mesh  for  the  sphere  used  to 
simulate  the  layered  immobile  region  was  manually  generated  by  assuming  the  region  from 
the  center  of  the  sphere  to  the  outside  region  would  consist  of  31  elements.  Each  node  of 
the  immobile  zone  finite  element  mesh  diffuses  through  its  boundaries  in  the  direction  of 
the  concentration  gradient.  Originally  the  plume  was  generated  using  a  non-sorbing  solute 
and  the  spheres  had  not  "charged  up"  with  contaminant.  The  first  time  step  of  the  initial 
pumping  simulation  "charged"  the  immobile  zone  by  setting  the  concentration  within  all 
spherical  mesh  elements  equal  to  the  initial  mobile  zone  concentration.  Thereafter,  the 
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immobile  zone  discharged  by  means  of  a  flux  toward  the  mobile  zone,  thereby  increasing 
the  mobile  zone  concentration  (causing  tailing  and  rebound)  and  decreasing  the  immobile 
zones  concentration  toward  an  equilibrium  state.  At  equilibrium  all  nodes  within  the 
immobile  zone  achieved  a  steady  state  concentration  and  the  flux  there  from  was 
necessarily  very  small  (eventually  approaching  zero).  The  finite-element  mesh  generated 
for  the  immobile  zone  was  fine  at  either  end  of  the  concentration  profile  (center  and 
exterior)  and  coarser  in  the  center  of  the  profile  where  gradients  were  not  as  large. 
Appendix  C  summarizes  the  finite-element  radii  for  the  layered  half-widths  for  the  Borden 
immobile  regions  studied. 

Boundary  Conditions 

Boundary  conditions  for  ground  water  flow  consisted  of  specified  head  boundaries 
on  all  sides  of  the  model,  as  depicted  in  Figure  4. 1  above.  Values  of  head  were  fixed  to 
maintain  a  constant  hydraulic  gradient  of  0.005  in  the  aquifer.  At  the  right  and  left  edges 
of  the  model,  as  indicated  in  Figure  4.1,  the  regional  ground  water  velocity  (0.079 
meters/day  [28.8  meters/year])  was  represented  by  a  pressure  value,  with  the  high  pressure 
on  the  left  and  a  low  pressure  value  on  the  right  edge  of  the  mesh.  No  flux  boundaries 
were  established  at  the  top  and  bottom  edges  of  the  model  respectively.  In  the  absence  of 
pumping,  these  boundary  conditions  simulated  spatially  constant  and  steady  ground  water 
flow  from  the  left  to  the  right  through  the  aquifer. 

Boundary  conditions  for  contaminant  transport  consisted  of  constant  zero 
concentrations  specified  at  all  model  boundaries.  All  water  entering  the  aquifer  from  the 
boundaries  was  assumed  contaminant  free.  Longitudinal  and  transverse  dispersivity 
parameters  were  chosen  in  an  iterative  manner,  vice  use  of  the  Borden  values,  to  ensure 
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that  the  gradient  did  not  migrate  near  the  previously  selected  mesh  boundaries.  However, 
the  longitudinal-to-transverse  dispersivity  ratio  (12:1)  was  within  the  range  of  values 
found  within  the  literature  as  being  representative  of  Borden  (2.5-20:1).  Selection  of  this 
operating  constraint  whereby  plume  migration  was  kept  within  the  mesh  boundaries 
negated  the  need  for  any  plume  contaminant  controls  during  pump-off  cycles  of  the 
remediation  design. 

Initial  Contaminant  Bums 

Steady  state  ground  water  flow  and  transient  solute  transport  were  simulated. 
Contaminant  was  simulated  by  injection  of  CTET  into  the  matrix  at  the  rate  of  5  m^/day, 
with  a  concentration  of  5.0  fig/1,  for  a  period  of  100-days.  While  the  initial  concentration 
appeared  to  be  small  and  was  equivalent  to  the  EPA  drinking  water  standard,  this  was 
chosen  to  allow  relatively  rapid  "clean-up"  of  the  hypothetical  aquifer.  Larger 
concentrations  should  take  longer  to  remediate,  but  would  display  the  same  general  effect 
of  the  lower  concentration  case,  but  over  a  longer  time  scale.  Upon  completion  of 
pumping  of  contaminant  into  the  matrix,  the  injection  well  was  turned  off  for  an  additional 
period  of  365-days.  Advective/dispersive  forces  in  conjunction  with  the  pressure  gradient 
(simulating  clean  water  flow  into  the  matrix)  continued  to  act  on  the  plume,  spreading  it 
through  the  simulated  mesh.  The  contaminant  plume  is  shown  in  Figure  4.2.  Solute 
concentrations  in  Figure  4.2  range  from  0  to  76  p.g/1.  The  elliptical  shape  and 
symmetry  of  the  plume  resulted  from  the  hypothetical  conditions  of  the  contaminant 
source,  homogeneous  flow  and  transport  in  the  simplified  aquifer  system.  This  plume 
served  as  a  tool  for  approximating  and  analyzing  remediation  strategies  incorporating  rate- 
limited  sorption  for  the  remainder  of  this  thesis. 
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CONCENTRATION  (g/1) 


CONCENTRATION  (g/l) 


n  Transverse  Dir.  (m) 
Down-Gradient  (m)  25  0 

Figure  4.2.  Initial  Contaminant  Plume  for  Hypothetical  1-D  Ground  Water  Model 


Aquifer  Remediation 
General 

The  remediation  system  consisted  of  a  single  ground  water  extraction  well,  located 
at  the  injection  point  previously  used  and  within  the  area  of  the  plume.  Pumping  was 
continuous  at  a  constant  rate  during  the  pump-on  cycle,  with  alternative  pump-off  (rest) 
cycles  until  the  aquifer  was  fully  remediated.  The  pumping  well  had  a  dual  role,  in  that  it 
not  only  extracted  the  contaminant  mass  from  the  aquifer,  but  also  functioned  as  a 
monitoring  well  for  cleanup  standard  achievement  by  mobile  zone  concentrations.  This 
dual  role  is  common  in  typical  field  cases.  Prior  to  performing  the  simulations,  it  was 
necessary  to  define  how  complete  remediation  of  the  aquifer  was  to  be  determined. 
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Pulsed  pumping  and  rest  period  criteria  also  required  definition.  After  these  criteria  were 
specified,  measures  of  optimality  or  efficiency  were  defined. 

Criteria  to  Define  Aquifer  Cleanup 

Water  quality  or  drinking  water  standards  are  sometimes  used  as  the  regulatory 
criteria  to  define  cleanup  of  an  aquifer  [Brogan,  1991:80].  A  rough  standard  for  drinking 
water  has  been  previously  defined  as  5.0  pg/1.  As  the  initial  concentration  was  arbitrarily 
low,  1%  of  this  value  (0.05  pg/1)  was  chosen  as  the  cleanup  standard  for  this  simulation. 

It  was  necessary  for  “rebound”  in  solute  concentrations  from  rate-limited  sorption  effects 
to  remain  below  this  criterion  after  pumping  was  stopped. 

LEA  Model 

For  unmodified  SUTRA  using  equilibrium  assumptions,  aquifer  cleanup  could  be 
achieved  when  the  solute  concentration  at  all  nodes  in  the  model  fell  below  0.05  p.g/1.  The 
output  files  from  SUTRA  facilitate  observation  of  the  concentration  at  all  nodes  within  the 
mesh.  Observation  was  made  of  the  entire  mesh  to  ensure  complete  capture  of  the  plume 
for  the  pumping  rates  selected  and  it  was  found  that  during  pumping  the  highest 
concentration  was  at  the  pumping  well  (due  to  inward  radial  flow).  When  the  cleanup 
criterion  was  met  and  no  rebound  was  observed  for  rate-limited  sorption,  the  remediation 
was  considered  complete,  pumping  was  stopped  and  the  simulation  ended.  Due  to  the 
equilibrium  assumption,  the  solute  concentration  would  continue  to  monotonically 
decrease  with  an  equilibrium  model  (due  to  lack  of  rebound  and  continued  advection  or 
extraction),  it  was  impossible  for  the  solute  concentration  to  rebound  following  the 
termination  of  pumping.  It  was  therefore  unnecessary  to  verify  that  the  equilibrium 
simulations  concentrations  remain  below  the  cleanup  criteria  after  cleanup  is  "achieved". 
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For  purpose  of  analysis  a  breakthrough  comparison  of  the  unmodified  SUTRA  was 
expected  to  show  earlier  breakthrough  and  achievement  of  the  cleanup  standard,  prior  to 
rate-limited  sorption  models,  as  well  as  a  monotone  decrease  with  pumping.  Continued 
pumping  beyond  the  cleanup  standard  was  moot,  in  this  case,  as  only  "clean"  water  would 
be  extracted. 

Non-Equilibrium  Model 

When  simulating  non-equilibrium  conditions,  local  equilibrium  between  the  sorbed 
and  mobile  phase  concentrations  was  not  assumed.  Although  mobile  (extracted)  phase 
concentrations  often  fell  below  the  cleanup  criteria,  contaminant  would  typically  reappear 
through  rebound  and  the  breakthrough  curve  would  exhibit  tailing  behavior  as  it 
asymptotically  sought  an  equilibrium  level.  Thus  for  both  pulsed  and  continuous  pumping 
situations,  cleanup  was  achieved  when  the  rebound  of  the  contaminant  no  longer  exceeded 
the  cleanup  standard.  The  standard  of  0.05  |J.g/l  (1%  of  drinking  water  standard)  was 
used  to  ensure  that  rebound  effects  from  desorbing  solute  did  not  exceed  the  primary 
drinking  water  standard. 

Comparing  Aquifer  Remediation  Strategies 

Two  primary  measures  were  used  for  comparison  of  different  remediation  designs: 
the  total  time  taken  to  achieve  cleanup  and  the  total  volume  of  contaminant  mass 
extracted  from  the  aquifer.  These  are  hereafter  referred  to  as  the  cleanup  time  and  the 
cleanup  mass  respectively.  Initially  a  pore  volume  (Vp)  analysis  was  conducted  to  assess 
the  normalized  volume  of  water  extracted  for  the  rate-limited  sorption  models  with  an 
LEA  model.  The  ultimate  measure  of  cleanup  effectiveness  would  include  an  economic 
analysis  of  remediation  costs.  During  the  entire  time  of  remediation,  for  instance,  there 
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are  costs  associated  with  operating  and  maintaining  the  ground  water  extraction,  treatment 
and  monitoring  systems,  in  addition  to  initial  capital  outlays  for  the  equipment  and 
associated  startup  costs.  Conceptually,  lower  pumping  rates  may  require  different 
equipment  and  thus  would  affect  the  overall  cost  of  the  system.  For  purpose  of  this 
analysis  all  capital  and  operating  costs  were  assumed  to  be  the  same.  Thus,  cost 
optimization  was  beyond  the  scope  of  this  thesis  and  efficiency  will  be  viewed  in  terms  of 
the  strategy  which  yields  the  greatest  mass  of  contaminant  during  the  pumping  period. 

This  strategy  was  chosen  for  the  following  reasons:  total  pumping  time  as  a  criterion 
would  require  attendant  economic  analysis  to  distinguish  amongst  different  pumping  rates 
(as  higher  pumping  rates  were  assumed  to  remediate  the  aquifer  faster  [not  necessarily 
more  efficiently]  and  would  always  be  superior  to  slower  pumping  rates);  total 
contaminant  removed  from  the  aquifer  which  was  the  ultimate  criteria  which  was  a 
consistent  factor  in  all  strategies  and  for  either  first-  or  second-order  non-equilibrium 
models. 


Non-Eauilibrium  Effects  on  Breakthromth  Curves 

First-  and  Second-Order  Pulsed  Pumping  Comparison 

Genual 

Pulsed-pumping  was  previously  defined  as  a  possible  improvement  on  continuous 
pumping.  Pulsed-pumping  is  attractive  in  that  it  would  allow  time  during  the  rest  cycle  for 
the  sorbed  solute  to  desorb  into  the  mobile  regions  for  later  extraction.  A  comparison  was 
made  between  the  mass  extraction  efficiency  and  breakthrough  arrival  times  for  the  two 
methods  of  modeling  physical,  rate-limited  sorption.  Rest  periods  for  the  "off-cycle"  of 
the  pumps  were  held  constant  and  were  reflective  of  the  inverse  ratio  of  the  diffusion 
coefficient  to  the  square  of  the  immobile  region  radius  [  l/(a^/De)  ].  This  criterion  is 
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similar,  but  of  lesser  magnitude  than  that  postulated  by  Goltz  and  Roberts  [1986a:87]  as 
represented  by  the  equation  1/a  .  Use  of  the  a^/De  ratio  provided  a  simple  means  of 
defining  rest  periods  for  both  first-  and  second-order  methods.  This  ratio  is  reflective  of 
the  time  scale  required  to  diffuse  contaminant  out  of  the  immobile  region.  First-order  rate 
constants  were  calculated  using  Equation  (3.29)  and  comparative  rest  periods  for  the 
immobile  zones  are  shown  in  Table  4.2  for  both  cases. 

Pumping  periods  were  variable  and  the  cessation  of  pumping  standard  was  based 
on  observation  of  the  time  at  which  the  mobile  region  concentration  decreased  below  the 
cleanup  standard  previously  defined  (0.05(1  g/1).  These  were  arrived  at  in  an  iterative 
manner  through  an  observational  approach,  whereas,  rest  periods  were  fixed  throughout 
the  simulation. 

A  total  of  1500  simulations  were  conducted  using  first-  and  second-order  cases  in 
the  range  of  Dg/a^  values  outlined  in  Table  4.2  for  pumping  rates  of  100  to  600  m^/day. 
Examination  of  the  extreme  cases,  for  both  small  and  large  De/a^  ratios,  was 
accomplished  through  breakthrough  curve  analysis  and  integration  of  portions  of 
breakthrough  curves  which  corresponded  to  pumping  periods  with  concentrations  above 
the  cleanup  standard  to  find  the  cleanup  mass  extracted.  Appendix  D  contains  the 
remainder  of  the  intermediate  breakthrough  curves  which  were  not  examined 
quantitatively. 


Pore  Volume  Analysis 

The  total  pore  volume  (Vp)  within  the  plume  was  estimated  by  the  equation 
[Valocchi,  1986:1689): 
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(4.5) 


Vp  =  Jifpfcn 


This  provided  an  estimate  of  the  total  time  required  to  pump  the  mobile  zone  free  of 
contaminated  water  and  replace  it  with  clean  water.  Using  Equation  (4.5)  the  pore 
volume  for  the  contaminant  plume  under  study  was  1620.31  m^.  Table  4.3  summarizes 
the  pore  volumes  required  to  treat  the  simulated  "sandy"  aquifer  using  an  equilibrium 
model  under  continuous  pumping,  and  first-  and  second-order  diffusion  models  under 
pulsed-pumping  conditions. 


TABLE  4.2 

COMPARATIVE  REST  PERIODS  FOR  Dg/a2  RATIOS  AND 
FIRST-ORDER  RATE  CONSTANTS 


fie 

a 

Ee/a2 

SL 

Tre  $  t_(a-/Pe) 

Trest  (1/SO 

(m2/d) 

(m) 

(d-1) 

(d-1) 

(d) 

(d) 

1.2  x  10*5 

0.05 

0.005 

5.8  x  10-8 

208.16 

198.67 

2.25  x  10'5 

0.05 

0.009 

1.1  x  10-7 

111.29 

106.21 

1.2  x  lO-5 

0.03 

0.013 

1.6  x  10-7 

74.94 

71.52 

3.5  x  10-5 

0.05 

0.014 

1.7  x  10-7 

71.44 

68.18 

2.25  x  10-5 

0.03 

0.025 

3.0  x  10-7 

40.06 

38.24 

3.5  x  10-5 

0.03 

0.039 

4.7  x  10-7 

25.72 

24.54 

1.2  x  10-5 

0.01 

0.12 

1.5  x  10-6 

8.32 

7.95 

2.25  x  10-5 

0.01 

0.225 

2.7  x  10-6 

4.45 

4.25 

3.5  x  IQ  5 

0.01 

0.35 

4.2  x  10-6 

2.86 

2.73 
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TABLE  4.3 


COMPARISON  OF  EXTRACTED  PORE  VOLUME  TO  ACHIEVE  CLEANUP 


EWa2 

Unmod 

Ratio 

SHEBA 

0=100 

Q=4QQ 

Q=6QQ 

0=100 

0.005 

2.3 

2.3 

2.4 

3.5 

0.009 

2.3 

2.3 

2.4 

2.9 

0.013 

2.3 

2.3 

2.4 

2.8 

0.014 

2.3 

2.3 

2.4 

2.8 

0.025 

2.3 

2.3 

2.4 

2.7 

0.039 

2.3 

2.3 

2.4 

2.7 

0.12 

2.3 

2.3 

2.4 

2.7 

0.225 

2.3 

2.3 

2.4 

2.7 

0.35 

2.3 

2.3 

2.4 

2.7 

lSL 

2nd- 

Order 

Order 

0=400 

0=600 

Q=1QQ 

Qf4QQ 

0=600 

7.0 

8.6 

2.8 

2.9 

3.0 

5.0 

6.2 

2.6 

2.9 

3.0 

4.1 

5.0 

2.6 

2.9 

3.0 

4.0 

4.9 

2.6 

2.9 

3.0 

3.2 

3.7 

2.4 

2.8 

2.9 

2.9 

3.1 

2.4 

2.7 

2.9 

2.6 

2.7 

2.4 

2.4 

2.5 

2.6 

2.7 

2.4 

2.3 

2.4 

2.6 

2.6 

2.3 

2.3 

2.4 

We  can  make  the  following  observations  from  the  data  in  Table  4.3: 

a.  Both  first-  and  second-order  diffusion  models  required  pumping  of  pore 
volumes  which  approached  that  of  the  equilibrium  model  (SUTRA  unmodified)  at  higher 
De/a^  ratios. 

b.  The  second-order  diffusion  model  was  more  precise  in  matching  the  extraction 
volume  attained  by  the  equilibrium  model  at  the  highest  Dg/a^  value. 
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c.  The  first-order  diffusion  model  required  treatment  of  more  pore-volumes  of 
ground  water  to  achieve  the  same  cleanup  level  as  the  second-order  model. 

d.  The  equilibrium  model  could  model  breakthrough  behavior  for  all  pumping 
rates  examined  at  the  higher  Dg/a^  values. 


Spherical  Diffusion 

The  second-order  diffusion  model  was  governed  by  the  flux  produced  by  the 
immobile  zone.  Figure  4.3  illustrates  the  gradients  at  various  points  in  time  for  the 
immobile  zone.  The  left  side  of  the  graph  represents  the  center  of  the  immobile  zone, 
while  the  right  side  of  the  graph  represents  the  exterior  of  the  sphere,  which  is  in 
equilibrium  with  the  mobile  zone  concentration.  The  upper  portion  of  the  graph 
represents  the  highest  concentrations  within  the  immobile  region,  while  the  lower  regions 
of  the  graph  represent  an  area  where  the  model  will  seek  equilibrium  with  the  mobile 
zones. 


At  the  start  of  pumping  the  immobile  zone  concentration  was  at  equilibrium  with 
the  mobile  fluid  concentration.  This  constant  concentration  throughout  the  zone  is 
depicted  with  the  linear  concentration  profile  at  the  initial  time  step  (T  =  0  days)  on  Figure 
4.3.  Upon  completion  of  the  first  pumping  cycle  (T  =  41.6  days),  the  immobile  zone  had 
discharged  contaminant  slowly  and  the  outer  elements  of  the  immobile  zone  mesh  had 
decreased  in  concentration,  while  the  interior  elements  remained  at  higher  levels. 
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Figure  4.3.  Low  Pumping  Rate  Immobile  Region  Concentration  Profiles 

(Q  =  100  m3/day)  (Dg/a2  =  0.005)  (Time  [T]  is  in  units  of  days) 


The  boundary  condition  (Equation  [3.27])  that  required  the  concentration  at 
element  31  (outer  node  of  immobile  zone)  to  be  in  equilibrium  with  the  mobile  zone  was 
confirmed  by  the  asymptotically  decreasing  concentration  gradient  which  terminated  at  the 
concentration  of  the  mobile  zone.  The  pump  was  then  turned  off  and  the  concentration 
profile  at  the  end  of  a  rest  cycle  (T  =  250-days)  was  approaching  an  intermediate, 
equilibrium  point,  as  depicted  in  Figure  4.3  by  the  small  slope  in  this  profile.  This  transient 
occurred  because  all  elements  within  the  immobile  zone  were  diffusing  or  absorbing 
concentration  fluxes  (as  contaminant  mass  moved  from  an  area  of  higher  concentration  to 
one  of  lower  concentration)  in  an  effort  to  equilibrate  at  the  same  concentration.  The  net 
result  of  this  equilibrium  seeking  behavior  was  a  decreasing  flux  out  of  the  immobile  zone, 
as  the  concentration  gradient  decreased  due  to  a  rise  in  the  mobile  zone  concentration  with 
subsequent  changes  in  the  outer  nodes  of  the  immobile  zone  mesh. 
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The  pump  was  then  turned  on  again  (due  to  the  rebound  in  concentration  observed 


within  the  mobile  zone  [as  shown  in  Figure  4.4]).  However,  as  the  concentration  of  the 
rebound  on  completion  of  the  first  rest  cycle  was  smaller  than  the  initial  concentration  of 
the  plume  (at  start  of  remediation),  it  took  a  shorter  time  to  reduce  the  concentration  of 
the  mobile  zone  below  the  action  level.  This  reduction  is  illustrated  by  Figure  4.4  below, 
where  at  completion  of  the  second  pump  period  (T  =  253.6-days)  the  BTC  once  again  is 
at  the  cleanup  level. 

Finally,  after  turning  the  pump  off,  once  again,  and  allowing  the  spheres  to  diffuse 
slowly  out  to  the  mobile  zone  for  another  rest  period  (T  =  400-days),  the  immobile  zones 
had  virtually  reached  equilibrium  once  again.  This  rest  period  was  less  than  the  standard 
period  in  an  effort  to  limit  the  focus  of  study  to  a  cleanup  time  of  less  than  400-days  for 
second-order  solutions  studied.  The  flux  out  of  the  immobile  region  was  very  minor  and 
did  not  cause  significant  rebound  within  the  mobile  region.  At  this  point  the  aquifer  was 
considered  cleaned  up,  as  the  mobile  region  concentration  was  below  the  cleanup  standard 
and  the  immobile  region  concentration  had  diffused  sufficiently  enough  that  it  would 
contribute  a  negligible  amount  of  contaminant  (T  =  400-days  in  Figure  4.3)  to  the  mobile 
region  for  rebound. 
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LEA 


Figure  4.4.  First-  and  Second-Order,  Pulsed  Pumping  Comparison  to  LEA  Model 
(Q  =  100  m3/day)  (Dg/a2  =  0.005) 

A  logarithmic  scale  is  used  to  enable  depiction  of  complete  BTCs.  A  linear  scale 
would  provide  only  a  limited  view  of  the  BTC.  This  relationship  is  illustrated  by 
Figure  4.5: 
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Figure  4.5.  First-  and  Second-Order,  Pulsed  Pumping  Comparison  to  LEA  Model 
(Q  =  100  m3/day)  (Dg/a2  =  0.005) 
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Small  Pe/a2  Ratio  Comparison  At  A  Higher  Pumping  Rate 

In  the  previous  analysis  (Figures  4.4  and  4.5)  a  small  Dg/a^  ratio  (0.005)  was 
analyzed  for  a  low  pumping  rate.  The  immobile  zone  layered  half-width  was  large  (0.05 
m)  and  the  effective  diffusive  rate  coefficient  (De)  [1.2  x  10' 5  m^/day]  was  small.  Now,  a 
larger  pumping  rate  (400  m^/day)  is  analyzed.  Figure  4.6  illustrates  these  points: 


Analysis  of  Figure  4.6  indicated  the  following: 

a.  The  first-  and  second-order  models  were  retarded  with  respect  to  the  LEA 
model  (represented  by  SUTRA  unmodified). 

b.  Both  first-  and  second-order  models  displayed  rebound  effects  after  the  pumps 
were  turned  off.  This  rebound  effect  was  due  to  diffusion  of  contaminant  from  the 
immobile  zone  into  the  mobile  fluid  by  means  of  concentration  gradient  fluxes. 
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c.  The  first-order  model  had  a  higher  rebound  level  than  the  second-order  model 
at  this  pumping  rate.  As  can  be  seen  in  Figure  4.4,  the  reverse  was  true  at  the  lower 
pumping  rate.  At  a  higher  pumping  rate  the  "shutoff  concentration"  in  the  mobile  zone  is 
reached  faster  than  with  the  lower  pumping  rate.  Therefore,  there  has  not  been  as  much 
time  for  diffusion  from  the  immobile  zone  during  the  pumping  cycle  and  thus  more  mass  is 
remaining  within  the  immobile  zone  in  both  models  for  the  higher  pumping  rate. 
Accordingly,  higher  rebound  thresholds  are  achieved  for  both  models  at  the  higher 
pumping  rate.  Apparently,  at  the  Q  =  lOOm^/day  pumping  rate,  during  the  first 
pumping  cycle  (T  =  41.6  days  in  Figure  4.4),  more  mass  moved  out  of  the  first  order 
immobile  zones  than  the  second-order  models  prediction  for  the  same  case.  With  the  Q 
=  400  m^/day  pumping  rate,  during  the  first  pumping  cycle  (T  =  22  days  in  Figure  4.6 
for  the  first-order  case)  the  reverse  is  true. 

d.  The  first-order  case  appeared  to  reach  equilibrium  faster  than  the  second-order 
case,  as  illustrated  by  the  sharp  flattening  in  its  concentration  profile.  The  second-order 
model  slowly  reached  its  equilibrium  point.  This  result  is  intuitive  as  the  immobile  zone 
nodes  seek  equilibrium  within  themselves  and  with  the  concentration  of  the  mobile  water 
surrounding  them,  acting  as  sinks  and  reservoirs  for  concentration  amongst  the  immobile 
zone  and  slow  the  rise  in  concentration  within  the  mobile  zone. 

e.  Departure  from  equilibrium  in  both  cases  was  rapid  and  evident  by  the  time  one 
pore  volume  of  ground  water  (4.05-days)  had  been  extracted.  This  result  was  also 
observed  by  Brogan  [1991:39]. 
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In  Figure  4.6  the  LEA  case  (SUTRA  unmodified)  appeared  to  achieve  the  same 
concentrations,  at  the  same  time,  as  the  second-order  case.  However  this  apparent 
convergence  was  due  to  the  graphical  resolution.  Figure  4.7  illustrates  this  point. 


Time  (Day) 

Figure  4.7.  SUTRA  Unmodified  Comparison  to  Second-Order  Diffusion  Model 
(Q  =  400  m3/day)  (Dg/a2  =  0.005) 


The  effect  of  an  even  higher  pumping  rate  (600  m^/day)  was  to  reduce  the  cleanup 
tifne  it  would  take  to  achieve  breakthrough,  yet  greater  pore  volumes  would  be  extracted 
to  achieve  cleanup  to  the  same  standard  within  the  aquifer.  Figure  4.8  illustrates  the 
breakthrough  curves  for  first  and  second  order  pulsed-pumping. 
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Figure  4.8.  First-  and  Second-Order,  Pulsed  Pumping  Comparison  to  LEA  Model 
(Q  =  600  m3/day)  (Dg/a2  =  0.005) 


Analysis  of  Figure  4.8  was  as  follows: 

a.  The  first-order  models  breakthrough  was  delayed  with  respect  to  the  second- 
order  model  (as  in  the  Q  =  400  m^/day  case)  due  to  the  relatively  constant  AC  flux 
inherent  in  the  first-order  models  algorithm. 

b.  The  second-order  model  achieved  breakthrough  earlier  due  to  the  limited 
diffusion  which  occurred  from  the  immobile  regions  during  the  time  the  pumping  was 
occurring.  This  result  is  expected  since  the  advective  time  scale  is  much  shorter  than  the 
diffusion  time  scale  for  this  case.  Furthermore,  physically  the  behavior  is  realistic  because 
higher  pumping  rates  flush  the  mobile  zone  clear  of  contaminant  with  clean  water  quickly 
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(faster  mobile  zone  velocity)  before  the  spherical  pore  elements  discharged  into  the  mobile 
zone  (slow  diffusive  velocity).  The  net  effect  within  the  immobile  regions  elements  was  to 
induce  a  sharper  concentration  gradient  in  concentration  within  the  outer  elements  of  the 
immobile  region  mesh  as  nodes  discharged  a  flux  to  maintain  the  boundary  conditions 
(Equations  [3.27]  and  [3.28]). 

c.  The  first-order  model  discharged  more  contaminant  from  the  immobile  zones  in 
retarding  the  breakthrough  curve  during  pumping,  yet  predicted  a  larger  rebound  after  the 
pumping  was  terminated.  This  illustrated  the  counter-intuitive  nature  of  the  first-order 
model  in  this  case. 

d.  First-order  diffusion  displayed  a  consistently  higher  rebound  than  second-order 
diffusion  at  the  higher  pumping  rate.  This  larger  rebound  would  result  in  more  pore 
volumes  removed  (by  the  first-order  model)  to  achieve  the  cleanup  standard  than  the  same 
second-order  case  (as  confirmed  by  Table  4.3). 

e.  First-order  models  reached  equilibrium  quicker  than  second-order  diffusion 
models  did. 

f.  After  one  pore  volume  (2.7-days)  had  been  removed,  first-order  models 
appeared  to  diverge  from  the  second-order  solution. 

g.  The  second  pulsed-pump  period  for  the  seccnd-order  model  was  less  than  the 
corresponding  first-order  model  pumping  period. 
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The  effect  of  the  larger  pumping  rate  within  the  immobile  region  was  to  induce  a 
sharper  gradient  in  the  concentration  profile  for  the  outer  elements  of  the  immobile  region 
mesh  as  nodes  discharged  a  flux  to  maintain  the  equilibrium  boundary  condition.  This 
gradient  was  illustrated  in  Figure  4.9  by  the  T  =  6.7  day  profile.  This  represented  the 
immobile  zone  concentration  upon  completion  of  pumping.  In  that  case,  the  inner 
immobile  region  elements  had  previously  maintained  concentration  at  higher  levels  and  the 
lower,  intermediate  equilibrium  point  reached  was  higher  than  that  in  the  100  m^/day 
pumping  case  illustrated  by  Figure  4.3.  When  the  pump  was  turned  on  again  a  greater 
pore  volume  extraction  (0.48)  was  required  to  reach  the  cleanup  standard  than  that  of  the 
lower  pumping-rate  case  (0.24).  The  net  effect  as  illustrated  in  Table  4.3  was  that  higher 
pumping  rates  must  extract  higher  pore  volumes  of  water  to  achieve  the  same  cleanup 
level  at  low  Dg/a^  values  as  lower  pumping  rates. 

Appendices  E  and  F  tabulate  the  pumping  strategies  required  to  achieve  cleanup  by 
both  the  first-  and  second-order  methods  for  all  pumping  rates  and  Dg/a^  ratios  examined. 
These  simulations  confirmed  that  after  a  second  pulsed  pump  period,  the  second-order 
model  did  not  display  any  significant  rebound  effects  for  the  low  and  high-pumping  rates 
previously  examined,  whereas  the  first-order  model  continued  to  require  successive 
pumping  and  rest  periods  to  achieve  the  cleanup  standard.  This  resulted  in  markedly 
different  cleanup  times  for  the  two  non-equilibrium  models  examined.  Table  4.4  illustrates 
the  difference  in  cleanup  time  for  the  pumping  rates  used  in  this  case. 
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TABLE  4.4 


CLEANUP  TIME  FOR  LOW  Dg/a2  RATIO  (0.005) 


Pumping  Rate 
(m^/day) 
100 
400 
600 


Cleanup  Time 
First-Order  Model 
(days) 

473 

643 

648 


Cleanup  Time 
Second-Order  Model 
(days) 

254 

220 

216 


Figure  4.9.  High  Pumping  Rate  Immobile  Region  Concentration  Profiles 

(Q  =  600  m^/day)  (Dg/a2  =  0.005)  (Time  [T]  is  in  units  of  days) 


Large  D^JaZ  Ratio  Case  Comparison 

In  the  case  of  a  large  Dg/a2  ratio  (0.35)  small  layered  half-widths  (0.01  m)  were 
assumed  to  have  a  large  effective  diffusion  coefficient  (De)  (3.5  x  10*5  m2/day). 
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Figures  4.10, 4.12  and  4.13  illustrate  the  breakthrough  curves  for  the  three  pumping  rates 
examined.  The  following  observations  were  made  by  analysis  of  the  breakthrough  curves 
and  Appendices  E  and  F  for  the  high  Dg/a^  ratio: 

a.  The  breakthrough  curves  for  all  the  cases  exhibited  little  difference  between  the 
first-  and  second-order  models  and  both  non-equilibrium  models  compared  closely  with 
the  LEA  model  until  the  initial  breakthrough  concentration  was  achieved.  Thereafter,  the 
non-equilibrium  models  exhibited  tailing,  while  the  LEA  model  continued  to  predict 
decreased  concentrations.  These  similarities  occurred  when  case  Peclet  numbers  for  the 
various  pumping  rates  were  the  largest  of  those  studied  and  were  approaching  the  limit  of 
the  non-equilibrium  region  where  molecular  effects  have  influence  on  dispersion  within  the 
mobile  zone.  This  result  illustrated  that  Peclet  number  may  not  be  as  valid  a  measurement 
of  non-equilibrium  effects  as  D^a2  . 

b.  The  first-order  model  predicted  minor  rebound  in  this  case,  which  was  not 
predicted  by  the  second-order  model  (as  found  in  Appendices  E  and  F).  The  negligible 
rebound,  due  to  rapid  dissipation  of  the  immobile  zone  solute,  was  confirmed  when  the 
immobile  region  concentration  profile  was  examined.  In  this  case,  the  large  diffusion 
coefficient  barely  retarded  movement  of  the  contaminant  out  of  the  immobile  region  into 
the  mobile  zone.  Thus  pumping  of  the  mobile  zone  served  to  remove  contaminant  equally 
as  well  from  the  immobile  region  due  to  the  smaller  diffusive  time  scale  (2.86-days  from 
Table  4.2).  At  the  end  of  the  first  pumping  period  (T  =  37.7  days  in  Figure  4. 1 1  and  T  = 
6.5  day  in  Figure  4.14)  the  contaminant  in  the  immobile  region  had  virtually  all  moved  into 
the  mobile  zone  and  there  was  almost  no  flux  of  contaminant  out  of  the  immobile  region 
to  facilitate  rebound,  due  to  depletion  of  mass  in  the  immobile  region  during  pumping. 
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c.  LEA  modeling  could  adequately  model  breakthrough  and  first-  or  second-order 
models  do  not  lend  greater  accuracy  in  results  for  the  computational  effort  required. 

d.  Pulsed  pumping  in  this  case  would  not  promote  greater  efficiency.  Continuous 
pumping  (at  all  pumping  rates)  for  this  Dg/a2  ratio,  removed  approximately  the  same 
number  of  pore  volumes  (from  Table  4.3)  of  ground  water  to  achieve  cleanup  of  the 
aquifer  as  pulsed-pumping.  Use  of  another  optimization  factor  (other  than  pore  volumes 
removed)  such  as  total  pumping  time  or  economic  factors  would  be  more  appropriate  in 
this  case.  Table  4.5  tabulates  the  cleanup  time  required  by  continuous  pumping  to 
remediate  the  aquifer  under  these  conditions. 


Figure  4.10.  First-  and  Second-Order,  Pulsed  Pumping  Comparison  to  LEA  Model 
(Q  =  100  m^/day)  (Dg/a2  =  0.35) 
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Immobile  Zone  Node  (from  ceiter) 


sure  4.1 1.  Low  Pumping  Rate  Immobile  Region  Concentration  Profiles 

(Q  =  100  m3/day)  (Dg/a2  =  0.35)  (Time  [ T]  is  in  units  of  days) 


Time  (Day) 

Figure  4.12.  First-  and  Second-Order,  Pulsed  Pumping  Comparison  to  LEA  Model 
(Q  =  400  m3/day)  (Dg/a2  =  0.35) 
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.  First-  and  Second-Order,  Pulsed  Pumping  Comparison  to  LEA  Model 
(Q  =  600  m^/day)  (Dg/a^  =  0.35) 


Inmobile  Zone  Node  (from  cater) 


Figure  4.14.  High  Pumping  Rate  Immobile  Region  Concentration  Profiles 
(Q  =  600  m-fyday)  (Dg/a^  =  0.35)  (Time  [T]  is  units  of  days) 


4-32 


TABLE  4.5 


CLEANUP  TIME  FOR  HIGH  Dg/a2  RATIO  (0.005) 


Pumping  Rate 
(m^/day) 


Cleanup  Time  Cleanup  Time 

First-Order  Model  Second-Order  Model 

(days)  (days) 


100  46.1  37.7 

400  13.5  9.5 

600  10.0  6.5 


Summary 

Up  to  this  point  we  have  been  solely  concerned  with  cleanup  time  and 
breakthrough  curve  analysis.  A  summary  of  what  has  been  discussed  so  far  is: 

a.  For  high  Dg/a2  ratios,  LEA  models  may  provide  an  efficient  means  of 
estimating  cleanup  time,  while  non-equilibrium  models  may  more  effectively  predict  tailing 
than  the  LEA  model. 

b.  If  the  objective  was  rapid  cleanup  and  an  LEA  model  applied,  then  higher 
pumping  rates  with  continuous  pumping  may  prove  to  be  the  most  effective  strategy. 

c.  Pulsed-pumping  at  lower  pumping  rates  resulted  in  a  longer  cleanup  time  for  the 
aquifer,  but  required  a  reduced  pore  volume  to  be  extracted  from  the  aquifer,  when  low 
Dg/a2  ratios  applied. 
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d.  First-order  models  may  over-  and  under-predict  rebound  within  the  aquifer, 
require  treatment  of  more  pore  volumes  and  predict  longer  cleanup  times  than  second- 
order  models  for  low  Dg/a^  ratios. 


EumBin&StcatceyXgmparisan 

We  now  examine  various  strategies  for  a  single  pumping  situation  using  the 
second-order  model  in  an  effort  to  determine  optimality  for  cleanup  mass  and  cleanup 
time.  This  analysis  examined  a  continuous  pumping  case,  a  pulsed-pumping  case,  and  a 
continuous,  effective  pumping  rate  case  for  a  single  Dg/a^  ratio  (0.005). 


A  continuous,  effective  pumping  rate  (Qeff)  can  be  hypothesized  once  the  pulsed- 
pumping  rates  are  determined  or  arrived  at  iteratively.  Brogan  [1991:137]  defined  the 
effective  pumping  rate  as  the  cleanup  volume  divided  by  the  cleanup  time: 


cu  = 


Q.*£t. 


pump 


y  t  +y  t 

^  pump  ^  re 


(4.6) 


The  continuous,  effective  pumping  rate  must  provide  for  containment  of  the  contaminant 
plume  to  be  considered  truly  effective.  Another  factor  in  its  formulation  was  that  the 
pulsed-pumping  rate  must  have  been  determined  already  to  facilitate  this  calculation. 
Finally,  the  continuous,  effective  pumping  rate  ensures  that  the  same  cleanup  volume  is 
extracted  from  the  aquifer  as  pulsed-pumping  methods.  However,  it  may  extract  less 
cleanup  mass  than  the  pulsed  method  (higher  rate).  Fickian  Diffusion  will  be  used  in  the 
following  calculations  as  was  discussed  in  Chapter  II  and  shown  in  this  Chapter.  First- 
order  methods  may  over  predict  mass  in  the  immobile  zones  at  high  pumping  rates  and 
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may  not  accurately  simulate  diffusion  in  high  flow  situations  or  in  cases  with  small 
sampling  scale  relative  to  the  diffusion  rate.  Finally,  second-order  (Fickian)  diffusion 
methods  are  more  firmly  founded  on  a  physical  basis.  Figure  4.15  illustrates  the  pumping 
strategies  that  will  be  compared  for  cleanup  efficiency  : 


Time  (D&y) 

Figure  4.15.  Pumping  Strategy  Comparison  (Dg/a^  =  0.005) 


Cleanup  mass  was  determined  by  composite  trapezoid  integration  [Burden  and 
Faires,  1993:187]  of  areas  under  the  concentration  curves  which  were  above  the  cleanup 
standard  for  the  various  strategies.  It  is  uneconomical  to  treat  water  which  was  pumped 
from  areas  with  concentrations  below  the  cleanup  standard  and  these  areas  were  not 
considered  to  contribute  to  the  overall  cleanup  mass  extracted  from  the  aquifer.  An 
additional  problem  with  this  method,  which  isn't  explored  is  what  becomes  of  the  treated 
water  (or  untreated  in  case  of  continuous  pumping  at  levels  below  the  cleanup  standard), 
as  some  authorities  may  not  allow  re-injection  of  the  pumped  water  into  the  aquifer  again. 
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The  total  cleanup  mass  is  arrived  at  by  use  of  the  following  equation: 


Md_=XpQ.^[C(»)+C(b)+22;f(xi)]  (4.7) 

i  2  j=i 

where  h  =  (b-a)/m ;  m  is  the  number  of  subintervals  between  the  pumping  times 
represented  by  b  =  time  ended  and  a  =  time  started;  xj  =  a  +  jh  for  each  j  =  0, 1 , ...,  m;  and 
n  is  the  pumping  period  number.  Based  on  the  injection  rate  of  5.0  m^/day  at  a 
concentration  of  5.0  |ig/l  for  100  days  it  was  estimated  that  5.0  x  10"2  kg  (assuming  all  is 
converted  to  dissolved  or  sorbed  solute)  of  contaminant  entered  the  ground  water. 

Table  4.6  outlines  the  mass  extracted  per  day  by  the  respective  pumping  strategies. 
While  the  continuous  pumping  strategy  achieves  breakthrough  at  the  same  time  as  the 
second-order  pulsed-pumping  case,  it  must  continue  to  extract  water  with  contaminant 
diluted  within  this  standard  (at  lower  concentrations)  to  effectively  clean  the  aquifer  up. 
The  continuous  pumping  strategy  required  extraction  of  a  larger  number  of  pore-volumes 
(22.3)  to  achieve  the  same  level  of  cleanup  as  that  achieved  by  the  pulsed-pumping 
strategy  (2.9).  This  is  the  result  of  the  slow  diffusion  of  contaminant  out  from  the 
immobile  zone  into  the  mobile  fluid.  In  essence  the  continuous  pumping  case  is  diluting 
clean  water  with  low  levels  of  contaminant  which  have  slowly  diffused  from  the  immobile 
zone  as  a  mass  flux. 
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TABLE  4.6 


MASS  EXTRACTION  PUMPING  STRATEGIES 
(D^a2  =  0.005) 


Mass  per  Dav 

Pore  Volumes 

Total  Cleanup 

Total 

Extracted 

Extracted 

Time 

PumDine  Time 

Strategy 

(kg/Day) 

(Days) 

(Days) 

Pulsed-Pumping 

1.3  x  10'3 

2.9 

219.5 

11.7 

Continuous 

2.3  x  lO'4 

22.3 

90.3 

90.3 

Pump 

Continuous, 

6.6  x  10'5 

2.9 

219.5 

219.5 

Effective 

Pumping 


The  following  conclusions  may  be  made  from  Tables  4.3  and  4.6;  and  Figure  4.15: 

a.  Pulsed-pumping  and  continuous,  effective  pumping  had  the  same  cleanup  time, 
but  different  cleanup  mass  extracted. 

b.  Continuous  pumping  had  a  significantly  longer  pumping  time  than  pulsed- 
pumping  as  opposed  to  cleanup  time,  where  the  reverse  is  true. 

c.  Terminating  the  continuous  pumping  case  when  the  cleanup  standard  was 
achieved  would  result  in  the  rebound  exhibited  by  the  pulsed  case  within  the  matrix,  thus 
prolonging  cleanup  time  beyond  that  predicted  by  the  continuous  case. 
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d.  Continuous  pumping  required  remediation  of  more  pore  volumes  (22.3)  of 
ground  water  than  either  continuous,  effective  pumping  or  pulsed-pumping  (2.9).  This 
would  have  a  major  effect  on  an  economic  analysis  because  of  the  large  volume  of  ground 
water  extracted  by  the  continuous  method  which  would  require  special  handling  and 
disposal,  if  current  practice  by  regulatory  authorities  continues  to  prohibit  re-injection  of 
pumped  water,  whether  treated  or  not.. 

e.  Continuous,  effective  pumping  extracted  less  cleanup  mass  than  either  pulsed- 
pumping  or  continuous  pumping  strategies. 

f.  Pulsed-pumping  was  the  most  efficient  method  in  extracting  mass  per  unit 
volume  pumped  or  per  day,  while  continuous  pumping  extracted  more  mass  (although  the 
majority  of  the  difference  between  this  method  and  pulsed  pumping  was  in  water  extracted 
at  concentrations  below  the  cleanup  standard). 

g.  Continuous  pumping  required  the  shortest  cleanup  time  of  the  three  strategies 
studied,  but  required  7.7  times  more  pumping  time  than  pulsed  pumping. 


4-38 


Yi  Conclusions  and  Recommendations 


This  concluding  chapter  draws  together  the  research  presented  in  the  previous 
chapters.  It  begins  with  an  overview  of  the  research  effort.  Next,  it  summarizes  the 
findings  of  this  study,  and  finally,  it  lists  recommendations  for  model  improvements  and 
for  follow-on  research. 

Overview 

The  focus  of  this  research  was  to  model  the  fate  and  transport  of  contaminated 
ground  water  under  conditions  of  rate-limited  sorption.  The  modification  of  the  SUTRA 
code  was  based  on  rugged  and  well  tested  code,  which  IRP  managers  can  use  to  better 
estimate  the  cleanup  time,  mass  extracted,  plume  retardation,  and  rebound  or  tailing 
effects  when  considering  a  pump-and-treat  extraction  system  for  remediation.  The  main 
thrust  of  the  model’s  evaluation  was  to  study  its  performance  under  a  range  of  realistic 
parameters  and  to  examine  the  differences  in  cleanup  criteria  for  LEA,  first-  and  second- 
order  diffusion  models.  In  the  latter  case,  demonstration  of  the  variability  in  predicted 
results  among  the  various  models  was  a  critical  objective.  This  objective  was 
accomplished  through  numerous  simulations  of  a  pulsed-pumping  extraction  system 
using  the  LEA  and  first-order  models  and  contrasting  their  performance  with  respect  to 
the  physically  and  computationally  more  accurate  second-order  model. 

The  literature  review  clearly  revealed  that  research  during  the  past  two  decades 
had  contributed  to  a  greater  understanding  of  rate-limited  sorption  at  the  laboratory  and  in 
limited  field  scale  experiments.  This  progress  had  greatly  enhanced  the  study  and 
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modeling  of  transport  processes  and  impacts  on  treatment  technologies  for  contaminated 
aquifers.  A  large  number  of  models  had  been  postulated  and  marketed  over  the  years. 
Models  that  emulated  rate-limited  diffusion  equations  were  found  to  be  primarily 
analytical  in  nature  and  highly  specialized  for  the  system  studied.  Two  other  authors  had 
modeled  chemical  non-equilibrium  processes.  However,  no  software  packages  have  been 
marketed  to  aid  IRP  managers  and  field  workers  in  modeling  remediation  projects  which 
incorporate  this  behavior.  Pulsed-pumping  had  been  the  subject  of  speculative  literature, 
concerning  its  applicability,  with  some  authors  asserting  its  non-effectiveness  over 
conventional  methods  of  pump-and-treat  systems. 

This  research  formulated  a  modification  to  the  SUTRA  code  for  physical,  non¬ 
equilibrium  rate-limited  sorption.  The  code  was  modified  to  incorporate  both  first-order 
and  second-order  rate-limited  sorption  methods.  The  modified  code  was  validated  for 
both  sorption  models  with  an  analytical  code  and  through  incorporating  previously  tested 
modules  and  routines.  An  extensive  series  of  simulations  was  conducted  to  evaluate 
performance  of  a  pulsed-pumping  extraction  system  in  remediation  of  a  hypothetical 
aquifer,  using  realistic  parameters.  Simulations  modeled  sorption  by  assuming  either 
equilibrium  or  rate-limited  conditions  existed  with  the  rate-limitation  described  by  either 
a  first-order  equation  or  by  Fickian  diffusion  of  contaminant  through  a  spherical 
immobile  region.  Additionally,  these  simulations  evaluated  cleanup  performance  for 
three  traditional  approaches  to  pump-and-treat  remediation  and  assessed  the  applicability 
of  this  code  and  remediation  strategies. 
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Summary  of  Findings 


Simulations  were  used  to  demonstrate  the  capability  and  viability  of  the  model 
modification.  In  addition,  a  comparison  of  this  model’s  simulations  to  those  of  another 
model,  which  had  recently  been  validated  against  other  codes,  was  conducted  with 
favorable  results.  Model  simulations  that  were  conducted  revealed  the  following: 


a.  A  pulsed-pump  operation  may  allow  for  more  efficient  removal  of  contaminant 
mass  per  unit  volume,  while  pumping  for  less  time  than  conventional  continuous 
pumping  methods.  The  model  clearly  demonstrated  the  rebound  and  tailing  behavior  for 
either  first-  or  second-order  (Fickian)  diffusion  models  which  are  characteristic  of  the 
slow  desorption  processes. 


b.  The  pulsed-pump  simulations  showed  lower  pumping  rates  required  less  pore 
volumes  of  ground  water  required  removal  from  the  aquifer  when  compared  to  higher 
pulsed,  pumping  rate  cases. 

c.  First-order  model  predictions  of  rebound  concentrations  and  cleanup  time  for 
aquifers  differ  significantly  from  that  predicted  by  second-order  methods  for  low  De/a2 
ratios  for  all  pumping  rates  evaluated. 


d.  First-order  models  reached  equilibrium  faster  than  the  second-order  models 
and  consequently  required  shorter  rest-cycles  than  the  second-order  model  to  diffuse 
contaminant  into  the  mobile  region. 
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e.  At  small  Dg/a^  ratios,  higher  pumping  rates  cleaned  up  the  aquifer  faster  than 


low  pumping  rates. 

f.  LEA  methods  could  model  breakthrough  behavior  at  high  De/a^  ratios  for  all 
pumping  rates  simulated,  since  no  significant  rebound  occurred,  although  first-order 
methods  would  be  adequate  for  modeling  the  tailing  behavior. 


g.  Continuous,  effective  pumping  did  not  extract  as  much  contaminant  mass  as 
either  continuous  pumping  or  pulsed  pumping  methods  and  required  more  total  pumping 
time  than  either  method  to  achieve  the  cleanup  level. 

h.  Continuous  pumping  at  a  high  rate  extracted  more  contaminant  mass  and 
achieved  cleanup  quicker  than  either  pulsed  or  continuous,  effective  pumping  methods,  in 
terms  of  total  remediation  time. 

Recommendations 

As  discussed  in  chapters  III  and  IV,  the  code  was  not  optimized  for 
computational  speed,  which  ultimately  impacted  on  the  size  of  aquifer  and  contaminant 
concentrations  modeled.  Therefore,  the  primary  emphasis  for  future  research  should  be 
to  increase  the  computational  speed  of  the  code  and  to  evaluate  its  performance  at  field 
scale  on  an  actual  remediation  site.  More  specifically,  emphasis  should  be  focused  in  the 
following  areas: 

a.  Evaluate  performance  of  the  modified  code  in  emulating  higher  concentration, 
contaminant  movement  and  retardation  during  a  field  tracer  test.  This  would  provide  a 
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much  needed  assessment  on  the  viability  of  the  model  and  its  accuracy  for  use  in 
modeling  other  IRP  sites. 

b.  Extend  the  code  to  model  multi-species  transport  in  a  2-D  aquifer. 

c.  Conduct  a  thorough  evaluation  of  pulsed-pumping  strategies  versus 
continuous  pumping  methods  in  optimization  of  pump  placement  and  optimal  pumping 
rates  or  rest  periods  for  maximum  efficiency. 

d.  Evaluate  model  performance  in  simulating  contaminant  transport  in  a 
heterogeneous  aquifer  with  spatially  variable  immobile  zone  radii,  diffusion  coefficients, 
porosities  and  hydraulic  conductivities. 

e.  Incorporate  multi-process  non-equilibrium  (MFNE)  routines  with  second-order 
diffusion  equations  within  SUTRA  to  comprehensively  model  behavior  of  solutes  which 
exhibit  both  physical  and  chemical  non-equilibrium. 

f.  Evaluate  performance  of  the  model  in  predicting  transport  of  a  sorbing  solute 
within  fractured  media  for  a  nuclear  or  chemical  waste  site. 
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Input  Parameter 


Value 


Qw  =QW 
b  =  B 
aL  =  ai 
0  -  THETA 
0m  =  THETAM 
0im  =  THETAIM 
ps  =  RHO 
a  =  ALPHA 
a  =  A 
D0  =  DO 
De  =  DE 
C0  =  F1 
Qm  =  CA 

CP 

f  =  EF 
Kd  =  KD 
4>  =  PHI 
rp  =  XF 1 

At 

Q 

Pw 

p 


1000  nrVd 


1.2  m 
0.33 
0.2838 
.0462 

1810  kg/m3 
5.8212  x  10‘8  s'1 
.05  m 

8.5  x  10"10  s-1 
1.17  x  10"5m2/d 


-1.0  x  10 


0.86 
28  m 
0.1  day 
46 

1000  kg/m3 
4.4  x  10"^  ms^/kg 
1.0  x  10'3  kg/m^s^ 


Parameter 


Value 


Qw 

rw 

b 

Ps 

Pw 

«L 

aj 

e 

9m 

0im 

Dm 

Dim 

k 

K 

dh/dl 

v0 

P 

P 

4> 

vp 

De 

a 

Ax 

At 

rP 


100-600  m3/d 
0.02  m 
20  m 

1810  kg/m3 
1 000  kg/m3 
1.2  m 
0.1  m 
0.33 
0.2838 
0.0462 

3.75  x  10'5  m2/d 

5.39  x  10"5  m2/d 

-12  2 
9.5  x  10  mz 

6  x  10’ 3  m/s 

0.005 

0.079  m/d 

4.8  x  10’10m2/N 
1.0  x  10"3  kg/m2s2 
0.86 

1620.31  m3 
1.2- 3.5  x  10’5  m2/d 
0.01  -  0.05  m 
1  m 

0.01  - 1.0  d 
8.85  m 
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Appendix  C:  Immobile  Region  Finite  Element  Mesh 


Node 

a  =  0,05  m 

a  =  0.03  m 

a  =  0.01  m 

1 

0.000 

0.000 

0.000 

2 

0.005 

0.003 

0.0010 

3 

0.010 

0.006 

0.002 

4 

0.015 

0.009 

0.003 

5 

0.020 

0.012 

0.004 

6 

0.025 

0.015 

0.005 

7 

0.030 

0.018 

0.006 

8 

0.035 

0.021 

0.007 

9 

0.040 

0.024 

0.008 

10 

0.042 

0.0252 

0.0084 

11 

0.043 

0.0258 

0.0086 

12 

0.045 

0.027 

0.0090 

13 

0.046 

0.0276 

0.0092 

14 

0.0465 

0.0279 

0.0093 

15 

0.0470 

0.0282 

0.0094 

16 

0.0475 

0.0285 

0.0095 

17 

0.0480 

0.0288 

0.0096 

18 

0.0485 

0.0291 

0.0097 

19 

0.0490 

0.0294 

0.0098 

20 

0.04925 

0.02955 

0.00985 

21 

0.04950 

0.02970 

0.00990 

22 

0.04955 

0.02973 

0.00991 

23 

0.04960 

0.02976 

0.00992 

24 

0.04965 

0.02979 

0.00993 

25 

0.04970 

0.02982 

0.00994 

26 

0.04975 

0.02985 

0.00995 

27 

0.04980 

0.02988 

0.00996 

28 

0.04985 

0.02991 

0.00997 

29 

0.04990 

0.02994 

0.00998 

30 

0.04995 

0.02997 

0.00999 

31 

0.05 

0.03 

0.01 
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Appendix  Bl  Breakthrough  Curves 
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Figure  D.2.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  400  m3/day)  (Dg/a2  =  0.009) 
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Figure  D.3.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  600  m3/day)  (Dg/a2  =  0.009) 
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Figure  D.5. 


Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  400  m3/day)  (Dg/a2  =  0.013) 


Figure  D.6.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  -600  m3/day)  (Dg/a2  =  0.0 1 3) 
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Figure  D.7.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  100  m3/day)  (D^a2  =  0.014) 


Figure  D.8.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  400  m3/day)  (Dg/a2  =  0.0 14) 
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Figure  D.9.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  600  m3/day)  (Dg/a2  =  0.014) 
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Figure  D.10.  Comparison  of  First-  and  Second  Order  Diffusion  Models  to  LEA 
(Q  =  100  m3/day)  (Dg/a2  =  0.025) 
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Figure  D.  1 1 .  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  400  m3/day)  (Dg/a2  =  0.025) 


Figure  D.12.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  600  m3/day)  (Dg/a2  =  0.025) 
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Figure  D.13.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  100  m3/day)  (D^a2  =  0.039) 
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Figure  D.14.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  400  m3/day)  (Dg/a2  =  0.039) 
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Figure  D.  15.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  600  m3/day)  (Dg/a2  =  0.039) 


0  50  100  150  200  250  300  350  400  450 

Time  (Day) 


Figure  D.16.  Comparison  of  First-  and  Second- Order  Diffusion  Models  to  LEA 
(Q  =  100  m3/day)  (Dg/a2  =  0. 1 2) 
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Figure  D.20.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  400  m3/day)  (Dg/a2  =  0.225) 
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Figure  D.14.  Comparison  of  First-  and  Second-Order  Diffusion  Models  to  LEA 
(Q  =  600  m3/day)  (Dg/a2  =  0.225) 
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Appendix  E:  First-Order  Pulsed  Pumping  Strategies 


1 .  Pulsed  pumping  times  and  rest  (pump  off)  times  are  given  in  days  for  the  various 
Dg/a2  ratios  and  pumping  rates  studied  for  the  first-order  diffusion  model.  Cleanup  Time 
is  the  total  time  to  remediate  the  hypothetical  aquifer  to  a  level  where  rebound  no  longer 
exceeds  the  action  level. 


Dg/a2 

0.005 

0.005 

0.005 

0.009 

0.009 

0.009 

Qw 

100 

400 

600 

100 

400 

600 

Pump  On 

54.11 

25.84 

21.16 

46.85 

18.64 

15.35 

Pump  Off 

208.2 

208.2 

208.2 

111.3 

111.3 

111.3 

Pump  On 

2.13 

1.82 

1.5 

0.4 

1.24 

1.06 

Pump  Off 

208.2 

208.2 

208.2 

111.3 

111.3 

111.3 

Pump  On 

0.05 

0.53 

0.51 

0.02 

0.29 

0.31 

Pump  Off 

0 

208.2 

208.2 

0 

111.3 

111.3 

Pump  On 

0 

0.04 

0.07 

0 

0.02 

0.04 

Pump  Off 

0 

0 

0 

0 

0 

0 

Total  On 

56.29 

28.23 

23.24 

47.27 

20.19 

16.76 

Cleanup 

Time 

472.69 

652.83 

647.84 

269.87 

354.09 

350.66 

Dg/a2 

0.013 

0.013 

0.013 

0.014 

0.014 

0.014 

Qw 

100 

400 

600 

100 

400 

600 

Pump  On 

45.2 

15.45 

12.5 

45.08 

15.13 

12.21 

Pump  Off 

75 

75 

75 

71.4 

71.4 

71.4 

Pump  On 

0.1 

0.91 

0.83 

0.09 

0.88 

0.81 

Pump  Off 

75 

75 

75 

71.4 

71.4 

71.4 

Pump  On 

0.02 

0.17 

0.21 

0.02 

0.16 

0.2 

Pump  Off 

0 

75 

75 

0 

71.4 

71.4 

Pump  On 

0 

0.02 

0.03 

0 

0.02 

0.02 

Pump  Off 

0 

0 

0 

0 

0 

0 

Total  On 

45.32 

16.55 

13.57 

45.19 

16.19 

13.24 

Cleanup 

Time 

195.32 

241.55 

238.57 

187.99 

230.39 

227.44 

E-l 


De/a2 

0.025 

0.025 

0.025 

0.039 

0.039 

0.039 

Qw 

100 

400 

600 

100 

400 

600 

Pump  On 

44.16 

12.38 

9.41 

43.08 

11.43 

8.19 

Pump  Off 

40.06 

40.06 

40.06 

25.72 

25.72 

25.72 

Pump  On 

0.07 

0.46 

0.52 

0.07 

0.19 

0.29 

Pump  Off 

40.06 

40.06 

40.06 

25.72 

25.72 

25.72 

Pump  On 

0.03 

0.04 

0.09 

0.03 

0.02 

0.04 

Pump  Off 

0 

0 

40.06 

0 

0 

0 

Pump  On 

0 

0 

0.02 

0 

0 

0 

Pump  Off 

0 

0 

0 

0 

0 

0 

Total  On 

44.26 

12.88 

10.04 

43.18 

11.64 

8.52 

Cleanup 

Time 

124.38 

93.00 

130.22 

94.62 

63.08 

59.96 

IVa2 

0.12 

0.12 

0.12 

0.225 

0.225 

0.225 

Qw 

100 

400 

600 

100 

400 

600 

Pump  On 

43.33 

10.79 

7.28 

43.11 

10.67 

7.16 

Pump  Off 

8.33 

8.33 

8.33 

4.45 

4.45 

4.45 

Pump  On 

0.05 

0.04 

0.05 

0.03 

0.03 

0.03 

Pump  Off 

8.33 

8.33 

0 

4.45 

0 

0 

Pump  On 

0.02 

0.02 

0 

0.02 

0 

0 

Pump  Off 

8.33 

0 

0 

0 

0 

0 

Pump  On 

0.02 

0 

0 

0 

0 

0 

Pump  Off 

0 

0 

0 

0 

0 

0 

Total  On 

43.42 

11.85 

8.33 

43.16 

10.7 

7.19 

Cleanup 

Time 

68.41 

28.15 

16.66 

52.06 

15.15 

11.64 

De/a2 

0.35 

0.35 

0.35 

Qw 

100 

400 

600 

Pump  On 

43.2 

10.59 

7.1 

Pump  Off 

2.86 

2.86 

2.86 

Pump  On 

0.03 

0.02 

0.02 

Total  On 

43.23 

10.61 

7.12 

Cleanup 

Time 

46.09 

13.47 

9.98 

E-2 


Appendix  F:  Second-Order  Pulsed  Pumping  Strategies 


1.  Pulsed  Pumping  times  and  rest  (pump  off)  times  are  given  in  days  for  the  various 
D^a2  ratios  and  pumping  rates  studied  for  the  second-order  diffusion  model.  Cleanup 
time  is  the  total  time  to  remediate  the  hypothetical  aquifer  to  a  level  where  rebound  no 
longer  exceeds  the  action  level. 


Dg/a^ 

0.005 

0.005 

0.005 

0.009 

0.009 

0.009 

Qw 

100 

400 

600 

100 

400 

600 

Pump  On 

41.1 

9.9 

6.7 

39.5 

10 

6.8 

Pump  Off 

208 

208 

208 

111 

111 

111 

Pump  On 

4.5 

1.8 

1.3 

3 

1.7 

1.2 

Pump  Off 

0 

0 

0 

111 

0 

0 

Pump  On 

0 

0 

0 

0.2 

0 

0 

Pump  Off 

0 

0 

0 

0 

0 

0 

Total  On 

45.6 

11.7 

11 

42.7 

11.7 

8 

Cleanup 

Time 

253.6 

219.7 

219.0 

264.7 

122.7 

119.0 

De/a2 

0.013 

0.013 

0.013 

0.014 

0.014 

0.014 

Qw 

100 

400 

600 

100 

400 

600 

Pump  On 

41 

10.1 

6.8 

40.1 

10.1 

6.8 

Pump  Off 

75 

75 

75 

71.5 

71.5 

71.5 

Pump  On 

1.7 

1.5 

1.2 

1.9 

1.5 

1.2 

Pump  Off 

0 

0 

0 

0 

0 

0 

Pump  On 

0 

0 

0 

0 

0 

0 

Pump  Off 

0 

0 

0 

0 

0 

0 

Total  On 

42.7 

11.6 

8 

42 

11.6 

8 

Cleanup 

117.7 

86.6 

83 

113.5 

83.1 

79.5 

Time 


F-l 


Dg/a2= 

0.025 

0.025 

0.025 

0.039 

0.039 

0.039 

Qw 

100 

400 

600 

100 

400 

600 

Pump  On 

39 

10.2 

6.9 

38.6 

10.1 

6.9 

Pump  Off 

40 

40 

40 

26 

26 

26 

Pump  On 

0.5 

1.2 

1 

0.2 

0.9 

0.8 

Pump  Off 

0 

0 

0 

0 

0 

0 

Pump  On 

0 

0 

0 

0 

0 

0 

Pump  Off 

0 

0 

0 

0 

0 

0 

Total  On 

39.5 

11.4 

7.9 

38.8 

11 

7.7 

Cleanup 

Time 

79.5 

51.4 

47.9 

64.8 

37.0 

33.7 

De/a2 

0.12 

0.12 

0.12 

0.225 

0.225 

0.225 

Qw 

100 

400 

600 

100 

400 

600 

Pump  On 

38.4 

9.6 

6.5 

38.3 

9.5 

6.5 

Pump  Off 

0 

8.3 

8.3 

0 

0 

0 

Pump  On 

0 

0.2 

0.3 

0 

0 

0 

Pump  Off 

0 

0 

0 

0 

0 

0 

Pump  On 

0 

0 

0 

0 

0 

0 

Pump  Off 

0 

0 

0 

0 

0 

0 

Total  On 

38.4 

9.8 

6.8 

38.3 

9.5 

6.5 

Cleanup 

Time 

38.4 

18.1 

15.1 

38.3 

9.5 

6.5 

De/a2 

0.35 

0.35 

0.35 

Qw 

100 

400 

600 

Pump  On 

37.7 

9.5 

6.5 

Pump  Off 

0 

0 

0 

Total  On 

37.7 

9.5 

6.5 

Cleanup 

37.7 

9.5 

6.5 

Time 


F-2 
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Tills  research' postulates  and  demonstrates  incorporating  rate-limited  sorption  effects  in  the  USGS  SUTRA 


code  for  cleanup  of  a  hypothetical  sandy  aquifer  by  pump-and-treat  remediation  methods.  Contaminant 
transport  is  assumed  to  be  affected  by  advection,  dispersion,  and  rate-limited  sorption/desorption.  Sorption  is 
assumed  to  be  either  equilibrium  or  rate-limited,  with  the  rate- limitation  described  by  either  a  first-order  law, 
or  by  Fickian  diffusion  of  contaminant  through  a  spherical  immobile  pore  region.  Solutions  are  arrived  at  by 
split  operator  methods  for  the  transport  and  one-dimensional  Galerkin  solutions  for  the  solute  concentration 
equations.  The  resulting  model  is  tested  against  an  analytical  Laplace  transform  model  for  both  first-order  and 
Fickian  diffusion  methods  in  a  radial  pumping  simulation.  Model  simulations  are  used  to  evaluate 
equilibrium,  first-order  and  Fickian  diffusion  effects  for  pulsed  and  continuous  pumping  solutions  within  a 
hypothetical  sandy  aquifer.  These  show  that  equilibrium  methods  under-predicted  rebound  while  first-order 
methods  may  both  under-and  over-predict  rebound  within  the  matrix  for  certain  regions  and  may  be  equivalent 
to  Fickian  diffusion  in  equilibrium  regimes  for  cleanup  time  prediction.  Model  simulations  are  then  used  to 
show  the  efficiency  of  pulsed  pumping  methods  in  cleanup  mass  extraction  per  pumped  volume  for  a 
contaminated  aquifer  pump-and-treat  remediation  activity  versus  more  conventional,  continuous  pumping 
methods. 


Physical  Non-Equilibrium,  Rate-Limited  Sorption,  Fickian  Diffusion,  * _  153 

SUTRA,  Pump-and-Treat,  Pulsed  Pumping  ' p  '  :  'rv 


1*7*  --•■.•-IT-  -J  .SIHOkUON  Vi  8  StC.UHiT  i  CLmSSIHCA  nON  (  i  >  StOjKllV  UaSSiti  .AUuiM  I  2y.  U  ■  -  u  \  Ur-  .avnAf 

Unclassified  1  Unclassified  ‘  Unclassified  «  UL 


